On certain matrix algebras related to quasi-Toeplitz matrices

Let $A_\alpha$ be the semi-infinite tridiagonal matrix having subdiagonal and superdiagonal unit entries, $(A_\alpha)_{11}=\alpha$, where $\alpha\in\mathbb C$, and zero elsewhere. A basis $\{P_0,P_1,P_2,\ldots\}$ of the linear space $\mathcal P_\alpha$ spanned by the powers of $A_\alpha$ is determined, where $P_0=I$, $P_n=T_n+H_n$, $T_n$ is the symmetric Toeplitz matrix having ones in the $n$th super- and sub-diagonal, zeros elsewhere, and $H_n$ is the Hankel matrix with first row $[\theta\alpha^{n-2}, \theta\alpha^{n-3}, \ldots, \theta, \alpha, 0, \ldots]$, where $\theta=\alpha^2-1$. The set $\mathcal P_\alpha$ is an algebra, and for $\alpha\in\{-1,0,1\}$, $H_n$ has only one nonzero anti-diagonal. This fact is exploited to provide a better representation of symmetric quasi-Toeplitz matrices $\mathcal {QT}_S$, where, instead of representing a generic matrix $A\in\mathcal{QT}_S$ as $A=T+K$, where $T$ is Toeplitz and $K$ is compact, it is represented as $A=P+H$, where $P\in\mathcal P_\alpha$ and $H$ is compact. It is shown experimentally that the matrix arithmetic obtained this way is much more effective than that implemented in the CQT-Toolbox of Numer.~Algo. 81(2):741--769, 2019.


Introduction
In this paper, we analyze structural and computational properties of the linear space Pα spanned by the powers of the semi-infinite tridiagonal matrix Aα = (a i,j ) i,j∈Z + , where a i,i+1 = a i+1,i = 1, a 1,1 = α and a i,j = 0 elsewhere, i.e., More specifically, let us denote a(z) = i∈Z a i z i the formal Laurent power series having coefficients a i ∈ C and let T (a) be the associated Toeplitz matrix defined by t i,j = a j−i .Recall that T (a) is symmetric if and only if a i = a −i for i ∈ Z + .We show that a semi-infinite matrix A is in Pα if and only if there exists a vector a = (a i ) i≥0 such that, for a(z) = a 0 + ∞ i=1 a i (z + z −1 ), one has A = T (a) + Hα(a), where Hα = (h i,j ), h i,j = η i+j−1 is the Hankel matrix whose first column η = (η i ) is such that , where v 1 = 1, v 2 = α, v 3+i = θα i , i = 0, 1, . .., and θ = α 2 − 1.Here, we assume conditions on α and a i under which the components η i have finite values.Observe that the above matrix is a Toeplitz matrix in Hessenberg form.Moreover, for α = 0 the matrix is tridiagonal skew-symmetric, while, if α = ±1 the matrix is bidiagonal and lower triangular.Moreover we show that if |α| ≤ 1 and the coefficients a i are such that ∞ i=1 |a i | < ∞, then the matrices in Pα have 2-norm and infinity norm bounded from above, so that Pα is a Banach algebra if endowed with • 2 or with • ∞.
If |α| > 1 and the coefficients a i are restricted to values such that the function a(z) = a 0 + ∞ i=1 a i (z i + z −i ) is analytic in the annulus A(r) = {z ∈ C | r −1 < |z| < r} for some r > |α|, then Pα is still a Banach algebra.
The above results are rephrased in the case where the matrix Aα is truncated to size n by setting an,n = β for β ∈ C. Let us denote by A α,β the matrix obtained this way so that In this case, the span P α,β of the powers A i α,β , i = 0, 1, . . ., n − 1, is a matrix algebra formed by matrices of the kind T (a) + Hα(a) + JH β (a)J, where J is the flip matrix having ones on the antidiagonal and zeros elsewhere.This algebra, for α, β ∈ {−1, 0, 1} includes the algebras investigated in [2,9,10,11,15].
As an application of these results, we provide a better representation of quasi-Toeplitz matrices having a symmetric Toeplitz part.We recall that a semi-infinite matrix A = (a i,j ) i,j∈Z + is quasi-Toeplitz (QT) if it can be written as A = T (a)+E A where T (a) is the Toeplitz part and E A is a compact linear operator from ℓ 2 in ℓ 2 , with [3].Quasi-Toeplitz matrices are encountered in queueing and network models, like the Tandem Jackson Network [18], that can be described by random walks in the quarter plane [3,5].In [6] a matrix arithmetic for QT matrices has been designed and analyzed.This software, implemented as the Matlab toolbox CQT-toolbox, has been used in [7] to solve certain quadratic matrix equations with QT coefficients encountered in the analysis of QBD processes [4], [16].
Here we show that representing a QT matrix A = T (a) + E A , where T (a) is symmetric, in the form Pα(a) + K A , for K A compact, is more convenient from the theoretical and the computational point of view.
In fact, since the set Pα is an algebra, we find that f (Pα(a)) = Pα(f (a)), for any rational function f (z) defined in the set {a(z) | |z| = 1}.This fact allows to reduce the computation of the matrix function f (Pα(a)) to the composition of two functions.This property is not satisfied by the matrices T (a) since symmetric Toeplitz matrices do not form a matrix algebra.
From the computational point of view, we compared the representation of QT matrices given in the CQT-Toolbox, designed in [6], with the new representation based on the algebra Pα.As test problems, we considered the solution of a matrix equation with QT coefficients performed by means of fixed point iterations, and the computation of the square root of a real symmetric QT matrix.In both cases, with the new representation we obtain substantial acceleration in the CPU time with respect to the standard representation given in [6].
The paper is organized as follows.In the next Section 2, we briefly recall some preliminary notions needed in our analysis.In Section 3, we introduce the linear space spanned by the powers of Aα and analyze its structural properties.In Section 4 we discuss some boundedness issues and give conditions in order that Pα is a Banach algebra.The cases where α ∈ {−1, 0, 1} are presented in Section 5.In Section 6 we show an application of these results to provide a better representation of QT matrices.Section 7 reports the results of some numerical experiments and finally Section 8 draws the conclusions.

Preliminaries
Denote T = {z ∈ C | |z| = 1} the unit circle in the complex plane, Z the set of relative integers, +∞} be the Wiener algebra.Given a(z) = i∈Z a i z i ∈ W, let T (a) = (t i,j ) i,j∈Z + be the semi-infinite Toeplitz matrix with entries t i,j = a j−i ∈ C, for i, j ∈ Z + .
We say that A is a quasi-Toeplitz matrix (in short QT-matrix) if A = T (a)+E A , where E A represents a compact linear operator from ℓ 2 into ℓ 2 , with We refer to T (a) as the Toeplitz part of A and to E A as the compact part of A, and we write E A = κ(A).The set of QT matrices is denoted by QT .
Denote by W S ⊂ W the set of functions a(z) = i∈Z a i z i ∈ W having symmetric coefficients, i.e., such that a i = a −i , and by QT S the set of QT matrices of the kind A = T (a) + E A where a ∈ W S .
Given QT matrices A = T (a) and denote with the same symbol v the semi-infinite column vector v = (v 1 , v 2 , . ..)T .Define H(v) = (v i+j−1 ) i,j∈Z + the semi-infinite Hankel matrix associated with v(z).Observe that the first column of H(v) coincides with the vector v.
Given a(z) = i∈Z a i z i ∈ W, define the following power series in the Wiener class: a + (z) = ∞ i=1 a i z i and a − (z) = ∞ i=1 a −i z i .We recall the following classical result, see for instance [8].= T (a n ).
In the following, we restrict our attention to the subset QT S .

A family of matrix algebras
For a given α ∈ C, denote where e 1 is the first column of the semi-infinite identity matrix I, and blanks denote null entries.The goal of this section is to provide a structural analysis of the matrices belonging to the linear space spanned by the powers A i α , i = 0, 1, . . ., n for any n ∈ Z + , that is, In this regard we will provide a different basis of this space that allows one to represent each matrix in the space as the sum of a banded Toeplitz and a Hankel matrix where only the Hankel component depends on α.
More specifically, we will introduce a new set {P 0,α , P 1,α , . . ., Pn,α}, of symmetric QT matrices having a better representation than A i α , and prove that this set is a basis of the linear space Sn,α.
Set θ = α 2 − 1, and define the polynomials in the variable z Denote with the same symbol hn the infinite column vector formed by the coefficients of the term z i , for i = 1, . . ., n, in the polynomial hn(z), filled with zeros.Denote also the semi-infinite Hankel matrix associated with hn whose first column is the vector hn.Finally, define and observe that P 1,α = Aα.For notational simplicity, if the dependence on α is clear from the context, we will write Pn in place of Pn,α, and Hn in place of Hn,α.
In particular, for n = 4 we have h 4 = (θα 2 , θα, θ, α, 0, . ..)T and , where blank entries stand for zeros.We need to prove a few preliminary results to show that the set {P i,α : i = 0, 1, . . ., n} is a basis of the linear space Sn,α.The first result that we need is the following identity, that is a direct consequence of the binomial expansion of (z + z −1 ) n : if n is even. ( Another useful property concerns the compact part Kn := κ(A n α ) of the QT matrix A n α , i.e., Kn = A n α − T ((z + z −1 ) n ).
Lemma 1 The compact part Kn of A n α is a Hankel matrix.Moreover, its first column kn = Kne 1 is such that ) n ) and both addends are symmetric, then Kn is symmetric.We prove that Kn is Hankel by induction on n.Clearly, the property is true for n = 1 where Kn = αe 1 e T 1 .Assume by inductive assumption that Kn is Hankel.Then, since K n+1 = κ(AαA n α ), a direct inspection shows that where 1 H(g n ), we find that Removing the first row of K n+1 we are left with the submatrix of AαKn obtained by removing its first row.This coincides with the product    that is, the sum of Kn plus the matrix obtained by shifting up the rows of Kn by two places.Since both the addends are Hankel matrices by the inductive assumption, then their sum is Hankel as well.Now, since K n+1 is symmetric and its submatrix obtained removing the first row is Hankel, then the whole matrix K n+1 is Hankel.
In order to prove (6), it is sufficient to multiply both sides of ( 7) to the right by e 1 .This manipulation yields To complete the proof, it is sufficient to determine the value of σn : is the constant term of the Laurent polynomial g(z) n .In view of ( 5), we have γn = 0 if n is odd and γn = ( n n/2 ) if n is even.On the other hand, δn := e T 1 H(g n )e 1 is the coefficient of z +z −1 in the Laurent polynomial g(z) n , and applying once again (5) we get δn = ( n ⌊n/2⌋ ) if n is odd, while δn = 0 if n is even.This proves ( 6) with σn = αγn − δn that completes the proof.⊓ ⊔ Let hn(z) be the function defined in (2), consider the matrix Hn = H(hn) which is the compact part of Pn = T (z n + z −n ) + Hn, and recall that hn = Hne 1 is the first column of Hn.
The following result relates the vectors hn and the matrix Aα.

. Then we have
where we have set h 0 := e 1 .
Proof The expressions for Aαh 1 and Aαh 2 follow by a direct inspection.For n ≥ 3, consider the ith component in both sides of equation (8).For i = 1, we have (Aαhn) . This shows that ( 8) is satisfied in the first component.Concerning the case i > 1, we may equivalently rewrite the condition (Aαhn is the semi-infinite shift matrix, and get ( Now, we are ready to prove the main result of this section. Theorem 2 Let Kn = κ(A n α ) and Hn = H(hn), where hn(z) is the function defined in (2).For the vectors hn = Hne 1 and kn = Kne 1 we have

Proof
We prove this identity by induction on n.For n = 1 we have For the inductive step, assume that the identity holds for n and prove it for n + 1.
By Lemma 1 and by the inductive assumption, we have If n is odd, applying Lemma 2 yields Now, recall that h 0 = e 1 , and σn = −( so that we get which completes the inductive step for n odd.Now consider the case where n is even.Then applying Lemma 2 to (9) yields Since ), applying (10) yields that is the desired equation since, for n even,
Proof We separately show that the Toeplitz part and the compact part in both sides of (11) coincide.Concerning the Toeplitz part, it is sufficient to prove the equality of the symbols, that is, which coincides with equation ( 5).Concerning the compact part, since κ(A n ) and κ(P n−2i ) are Hankel matrices defined by their first column kn and h n−2i , respectively, we may equivalently rewrite the condition The latter vector equation follows directly from Theorem 2. ⊓ ⊔ The following result shows that equation (11) allows to express Pn,α as a linear combination of A i α , for i = 0, . . ., n.
Proof Clearly, equation (11) provides an explicit expression of A n α in terms of the matrices P i,α , i = 0, 1, . . ., n. Conversely, we may write (11) as By means of an inductive argument it follows that any matrix P j,α , for j = 0, 1, . . ., n, can be expressed as a linear combination of the matrices A i α for i = 0, 1, . . ., j.This completes the proof.⊓ ⊔ Given a symmetric Laurent polynomial a(z) = a 0 + n i=1 a i (z i + z −i ), denote Pα(a) ∈ Sn,α the matrix We say that Pα(a) is the QT matrix in Sn,α associated with the symbol a(z).As a synthesis of the results of this section we may conclude that That is, the symbol a(z) uniquely defines the Toeplitz part T (a) and the compact (Hankel) part Hα(a) of Pα(a).In particular, any matrix in Sn,α is the sum of a symmetric Toeplitz matrix and a Hankel matrix associated with the same symbol a(z).
One may wonder if the matrix Pα(a) can be defined in the case where a We provide an answer to this question in the next section.

Boundedness issues
A natural issue concerns the boundedness of in the case where a(z In order to answer this question, observe that from (2) we have It is interesting to point out that for |α| ≤ 1, the above quantity is bounded from above by 1 + 2|α| which is independent of n.
If |α| > 1, in general for a(z) ∈ W S it is not guaranteed that Pα(a) ∞ < ∞, in fact, according to (14), the norm Hn,α ∞ grows as O(|α| n ).However, assuming that a(z) ∈ W S is also analytic for z in the annulus for some r > |α| > 1, then it is possible to show the boundedness of Pα(a) ∞.In fact, from the analyticity of a(z) in A(r) it follows that the series ∞ n=0 |α| n |an| is convergent, so that ∞ n=1 anHn,α ∞ is bounded from above [12, Theorem 4.4c].A similar analysis can be carried out by replacing the infinity norm with the 2-norm.In fact, by definition, we have Hn,α 2 2 = ρ(H * n,α Hn,α), where ρ(•) denotes the spectral radius of a matrix and the superscript " * " denotes conjugate transposition.Since Hn,α is possibly nonzero only in the n × n leading principal submatrix Hn,α of Hn,α, we have Hn,α 2 2 = ρ( H * n,α Hn,α).Moreover we have 2 , where the matrix inequality is meant entry-wise.From the theory of nonnegative matrices [1], we know that if A and B are n × n matrices such that |B| ≤ A then ρ(B) ≤ ρ(A), thus we obtain where the latter inequality holds since ρ(A) ≤ A ∞. Thus, since Hn,α ∞ = Hn,α ∞, from ( 14) we get This inequality, together with the property T (a is analytic in the annulus A(r) defined in (15), for some r > |α| > 1.
Define the following set We may synthesize the analysis of this section by the following result.
and let Pα(a) be the matrix defined in (13).Assume that one of the following conditions is satisfied: |α| > 1 and a(z) is analytic in the annulus A(r) defined in (15), where r > |α|.As a consequence of the above properties we may say that if f (z) : C → C is a rational function defined in the set a(T), so that f (a(z)) ∈ W S , then f (Pα(a)) = P (f(a)) if |α| ≤ 1.Moreover, the same equality holds if |α| > 1 and f (a(z)) is analytic in A(r) for some r > |α|.
These properties reduce the matrix arithmetic in the class Pα to function arithmetic in W S , that can be easily implemented whatever is the value of α with |α| ≤ 1, by following the same techniques developed in [6] for QT arithmetic.The same arithmetic can be implemented for |α| > 1 in the case where the symbols are analytic in A(r).
From the algorithmic point of view, we restrict the attention to the case where ) are Laurent polynomials.In the case of Laurent series we need to truncate a(z) and b(z) to a finite degree.In this case, the operation of addition is almost immediate since for c(z) = a(z) + b(z) we have c i = a i + b i , i = 0, . . ., max(m, n) where we assume that a i = 0 and b i = 0 for i out of range.
The case where c(z) = a(z)b(z) is slightly more complicated.Here, for a given integer N we need a principal N -th root of 1, ω N = cos(2π/N)+i sin(2π/N), where i is the imaginary unit such that i 2 = −1.The Laurent polynomial c(z) is such that z m+n c(z) is a polynomial of degree at most 2(m + n).In order to compute the coefficients c i , it is sufficient to fix an integer N > 2(n + m), and compute the values a(ω i N ), i = 1, . . ., N with 2 FFTs, and then compute the N products The values of the coefficients c i are finally recovered by interpolating the pairs (ω i N , y i ) by means of IFFT.The overall cost of this procedure is O((m + n) log(m + n)) arithmetic operations.
A similar approach can be applied to compute Pα(a) −1 = Pα(1/a(z)) in the case where a(z) = 0 for |z| = 1.In fact, in this case the function 1/a(z) is a Laurent series that can be approximated by a Laurent polynomial c(z) of suitable degree.The coefficients and the degree of c(z), together with the ratio cond = max z∈T |a(z)|/ min z∈T |a(z)|, can be determined dynamically by means of Algorithm 1.
The value of "cond" provides an estimate of the numerical invertibility of a(z).Indeed, if cond=∞ the function is not invertible and if cond> ǫ −1 , where ǫ is the machine precision, the problem of inverting a(z) is numerically ill-conditioned.The cost of the above algorithm is O(N log N ) arithmetic operations, where N is the numerical degree of the Laurent series 1/a(z).

Some special cases
By choosing specific values for α we obtain very special subalgebras of the set of QT matrices having symmetric Toeplitz part.

The finite case
The analysis performed in Section 3 for semi-infinite matrices can be easily applied to m × m matrices.More specifically, from Theorem 2 we may deduce a structural characterization of the matrix algebra generated by the powers of the m×m matrix This algebra has been investigated in [9,10,11,15] in the cases where α, β ∈ {0, 1, −1}.
In fact, by adapting to the finite case the same argument used to prove Theorem 2 and its corollaries, we find that where i is the m × m leading principal submatrix of the semi-infinite Hankel matrix H i,α and Jm is the m × m permutation matrix with ones in the anti-diagonal.An example of the m × m matrix P (α,β) i is given for i = 4 and m = 8: For α = β = 0 we obtain the τ algebra analyzed in [2], where the matrices are simultaneously diagonalized by the discrete sine transform of kind 1.For α, β ∈ {−1, 0, 1} the corresponding matrices are simultaneously diagonalized by different types of discrete sine and discrete cosine transforms, see [15] for more details.These classes have been used to precondition positive definite Toeplitz systems.
For any value of α and β, the matrices in the algebras generated this way can be written as the sum of a Toeplitz and a Hankel matrix.

QT matrix representation
The implementation of QT matrix arithmetic given in the package CQT-Toolbox of [6] relies on the fact that, given two QT-matrices still compact.In fact, since the matrices are finitely representable, then E A and E B have finite rank so that also rank(E C ) is finite since rank(E C ) ≤ rank(E A ) + rank(E B ). Similarly, relying on Theorem 1, the product C = AB, is written as For this reason, a strategy of compression and approximate representation is introduced in [6] in order to reduce the growth of the rank when many arithmetic matrix operations are performed.This compression strategy is the most time consuming part of the software CQT-Toolbox [6].
Assuming that the decomposition H are available, where matrices U A , V A are formed by k A columns, U B , V B are formed by k B columns, and U H , V H are formed by k H columns, the compression strategy requires the computation of the numerical rank of the matrix In order to perform the compression, two matrices U C , V C having less than k A + k B + k H columns are computed such that E C is approximatively given by U C V T C , up to an error of the order of the machine precision.This computation, usually performed by means of SVD, is time consuming especially in the cases where the numerical size of the compact correction gets large.This happens, say, when the bandwidth of the Toeplitz part grows large in any algorithm using the package CQT-Toolbox.
Concerning matrix inversion, in [6] the inverse of A = T (a) where the inverse of I+T (a) −1 E A is computed by means of the Sherman-Woodbury-Morrison (SWM) formula, and the inverse of the Toeplitz matrix T (a) is computed relying either on the cylic reduction algorithm, or on FFT.For more details we refer the reader to [6].
For QT matrices A = T (a) + E A ∈ QT S , where a(z) ∈ W S , the results of Section 3 suggest us to represent A in the form where Pα(a) = T (a) + H(a) ∈ Pα, and In fact in this case we find that, for A = Pα(a) + K A , B = Pα(b) + K B , we have and c(z) = a(z)b(z).Observe that, with this representation, we simply have and K B = ÛB V T B .A similar advantage is encountered in the computation of A(a) −1 = (I + Pα(a −1 )K A ) −1 Pα(a −1 ).In fact, in this case, the computation of a(z) −1 , performed by means of Algorithm 1, is much easier than inverting T (a) and the major computational effort is just applying the SWM formula to the matrix I + Pα(a −1 )K A .
In general, if f (x) is a rational function defined on the set a(T), then f (A) = Pα(f (a)) + K f , where K f is compact.Moreover, in a straight-line program whose generic step is of the kind Sm = S i • S j where • ∈ {+, =, * , /} is any of the four operations, i, j < m and S ℓ = Pα(s ℓ )+K ℓ ∈ QT S , one finds that Sm = Pα(sm)+Km, with sm(z) = s i (z) • s j (z).Moreover, the compact correction Km is related to the compact corrections K i and K j of the operands and to the specific operation • used at the m-th step.
Instead, by using the representation Sm = T (sm) + Em given in [6], one finds that Sm = T (s i •s j )+Em, where this time Em is related not only to the corrections E i and E j of the operands, but also depends on the correction T (s i ) • T (s j ) − T (s i • s j ) that, in the case of multiplication or inversion, can have large rank.
In other words, one expects that the matrices to be compressed in the representation of [6] have rank larger than the rank of the matrices to be compress in the new representation through Pα.
The actual advantages of the new representation will be shown in the next section concerning the numerical experiments.

Numerical experiments
We say that A = Pα(a) + K A ∈ QT S is finitely representable if the symbol a(z) is a Laurent polynomial and the correction K A has a finite support, that is, its entries are zero outside its leading principal m × n submatrix for suitable finite values of m and n.
We have implemented the representation of a matrix in QT S in terms of the pair (Pα(a), K A ) together with the basic arithmetic operations.
In this section, we report the results of a preliminary experimentation that we have performed in order to compare the efficiency of the above representation with the standard representation of general QT matrices given in the CQT-Toolbox of [6], see https://github.com/numpi/cqt-toolbox.The experiments have been performed in Matlab v. 9.5.0.1033004 (R2018b) on a laptop with an Intel I3 processor and 16GB RAM.The functions and the scripts that execute the tests can be downloaded from https://numpi.dm.unipi.it/scientific-computing-libraries/sqt/ The first test concerns the solution of the quadratic matrix equation We may see how the new representation provides a speed-up in the CPU time that ranges roughly from 1.6 to 5.8.Moreover, comparing with the representation of X i as matrix in P 1 , we have speedups in the range [124,320].Observe also that the residual errors are close to the machine precision.
We may observe a slight discrepancy between the number of iterations required for the two different representations.This is due to the fact that the stop condition involves the norm of the compact correction that can be different in the two representations.It must be said that in the experiments involving the CQT-toolbox, cpu time, number of iterations, and errors can slightly change in different executions of the code.The reason is that the compression operation performed in the CQT-toolbox is randomized.
The second test concerns the computation of the square root of a matrix A performed by means of the algorithm [13, Equation (6.20)] designed in [14].We have considered the class of matrices given by Aγ = T (a), where T (a) is the symmetric Toeplitz matrix associated with the symbol a(z) = γ + 4(z + z −1 ) + 3(z 2 + z −2 ) + 2(z 3 + z −3 ) + (z 4 + z −4 )], for γ > 5.This matrix is positive definite for γ > 5, while if γ = 5 it is not invertible.Table 2 reports the CPU time in seconds, the number of coefficients of the symbol, the size and the rank of the correction, together with the residual error X 2 −A ∞, where X is the computed matrix, for values of γ = 5+δ with δ = 10 −1 , 10 −2 , 10 −3 .Between parentheses the CPU time required for computing the symbol is reported.This time corresponds to computing the square root of any matrix of the kind Pα(a).
Observe that also in this case, the speed up with respect to the QT representation ranges from 2.22 to 3.15.While comparing the CPU time for computing the symbol, i.e., for computing the square root of A = Pα(a), with the time required by the QT representation we get speed-ups in the range [78,233].
In all the experiments that we have performed, the residual errors obtained with the new representation are comparable with ones obtained with the CQT-toolbox.

Conclusions
A class, depending on a parameter α, of subalgebras of semi-infinite quasi-Toeplitz matrices having a symmetric symbol has been introduced.A suitable basis of these Table 2 Computation of the square root of the SQT matrix A = T (a), where a(z) = 5 + δ + 4(z +z −1 )+3(z −2 +z 2 )+2(z −3 +z 3 )+(z −4 +z 4 ), using the standard (QT ) and the new (QT S ) representation of A. Values of the CPU time in seconds, number of iterations, length of the symbol, size and rank of the compact correction, together with the residual error X 2 − A ∞, are reported.Between parentheses, in bold, the CPU time required for A = P 1 (a).
subalgebras, as linear spaces, has been determined and exploited to represent in a more efficient way quasi-Toeplitz matrices associated with a symmetric symbol.A preliminary experimentation has been performed concerning the computation of the matrix square root and the solution of a quadratic matrix equation.From our tests it turns out that the new representation is more effective with respect to the one given in the Matlab toolbox of [6].A structured analysis of these algebras has been carried out both for semi-infinite matrices and for finite matrices where the algebras depend on two parameters α and β.The latter class includes the algebras investigated in [2,9,10,11,15] obtained for α, β ∈ {−1, 0, 1}.To our knowledge, the case of general values of α and β has not been analyzed in the literature.

Theorem 1
For any a(z) ∈ W, the matrices H(a + ) and H(a − ) are compact operators.Moreover, for b(z) ∈ W it holds that T (a)T (b) = T (ab) − H(a − )H(b + ), where H(a − )H(b + ) is compact.The above result implies that for any QT matrices A = T (a) + E A , B = T (b) + E B it holds AB .= T (ab).As a consequence we have (T (a) + E A ) n .
and the Hankel part H(a − )H(b + ) disappears since it is included in Pα(c).This way, the compression of the correction requires the computation of the numerical rank of the matrix [ ÛA , Pα(a) ÛB ][Pα(b) VA + VB ( ÛT B VA ), VB ] T , where K A = ÛA V T A

Table 1
[7]negative matrices in QT S such that A + B + C is stochastic.In this case, the solution of interest is the minimal nonnegative solution G that exists and is a QT matrix, i.e., G = T (g) + E G[7].This solution can be computed by means of the classical Solution of the quadratic matrix equation (20) by means of Algorithm 1, 2, and 3 of (21) where the coefficients are QT matrices in P 1 , represented either as entries in P 1 (CPU time between parentheses and in bold), or in the form W + K for W ∈ P 0 and K compact (column QT S ), or in the form T + E, where T is Toeplitz and E is compact (column QT ).Besides the CPU time in seconds, the table reports the number of iterations, the number of coefficients in the symbol, the size and the rank of the correction, together with the infinity norm of the residual error AG 2 + BG + C − G ∞, where G is the computed solution.