Stability analysis of the matrix-vector product (via FFT) for a
Toeplitz-like matrix
PAOLA FAVATI1, ORNELLA MENCHI2
1IIT - CNR Via G. Moruzzi 1, 56124 Pisa, ITALY
2Dip. di Informatica, University of Pisa, Largo Pontecorvo 3, 56127 Pisa, ITALY
Abstract: -In this paper the numerical stability of the matrix-vector product, performed via FFT, is analyzed for
the case of the Toeplitz-like matrices. The error appears to depend on the magnitude of the generators of the
matrix. The numerical experimentation confirms the theoretical result.
Key-Words: -Stability analysis, Toeplitz-like matrix, Fast Fourier Transform
Received: May 4, 2021. Revised: January 11, 2022. Accepted: January 28, 2022. Published: February 25, 2022.
1Introduction
Let us consider a Fredholm integral equation of the
first kind
g(s) = ZK(s, t)x(t)dt, (1)
where the square integrable kernel K(s, t)and the
right-hand side g(s)are given functions and x(t)is
the unknown solution to be reconstructed. In many
applications g(s)consists of measured quantities.
Problems modeled by equation (1) are frequent both
in the one-dimensional context (for example in signal
processing and in the computation of inverse trans-
formations) and in the two-dimensional context (for
example in the image deconvolution problem where
K(s, t)represents an imaging system, x(t)and g(s)
represent a real object and its image, respectively).
By discretizing (1), a linear system
Ax=g(2)
is obtained, whose main features are the large dimen-
sion of Aand the distribution of its singular values
which often decay gradually to zero. Hence solving
(2) is an ill-posed problem. When some structure can
be considered for A, system (2) results to be solv-
able in practice also for large dimensions. In many
applications K(s, t)can be assumed to be invariant
with respect to translations and with a bounded sup-
port, so that matrix Aresults to have Toeplitz struc-
ture and a limited bandwidth. Toeplitz systems arise
frequently in linear algebra problems. Unfortunately,
when a simple operation like multiplication or inver-
sion or low rank modification is applied to a Toeplitz
matrix, the Toeplitz structure is lost. For example,
an important quantity as the Schur complement of a
leading principal submatrix of a Toeplitz matrix has
no longer a Toeplitz structure. For this reason we
consider here a more general structure, the class of
Toeplitz-like matrices, which is closed for the most
common operations applied in the numerical algo-
rithms. The Toeplitz-like structure is based on the
concept of displacement rank [1, 2, 3, 4] and has been
studied by many authors with applications in several
fields (see for example [5] for an extensive bibliog-
raphy). In the last decade several papers dealt with
fast and superfast solver (see for example [6, 7, 8],
application to the solution of fractional partial differ-
ential equations [9, 10, 11] and to the queueing prob-
lem [12]. Moreover a great interest was addressed to
the study of Toeplitz-like operators in infinite dimen-
sional spaces ([13, 14, 15])
The special low cost algorithms based on the fast
Fourier transform (FFT), which have been devised
to perform the matrix-vector products with Toeplitz
matrices, can be applied also to Toeplitz-like matrices.
FFT was first discussed by Cooley and Tukey in
1965 [16], although Gauss had already described the
critical factorization step as early as 1805. It is one
of the most important numerical algorithms and has
a wide range of applications. Its ubiquitous fortune
is principally due to a low computational cost: com-
puting the discrete Fourier transform of a sequence
of length naccording to the definition, takes O(n2)
arithmetical operations, while using FFT it takes only
O(nlog n)operations.
When finite-precision floating-point arithmetic is
used, FFT algorithms give results affected by error,
but this error is typically quite small, in fact most
FFT algorithms enjoy excellent numerical stability
(see [17, 18]). We are interested in investigating the
stability of the matrix-vector product based on FFT
for the case of Toeplitz-like matrices. The paper is so
organized: in Section 2 two simple programs describe
the use of FFT in the computation of the matrix-vector
product for triangular Toeplitz matrices, in Section 3 a
brief description of Toeplitz-like matrices is given and
the function which computes the matrix-vector prod-
uct for Toeplitz-like matrices is sketched. The analy-
sis of the stability occupies Section 4 giving an upper
bound of the error which depends on the magnitude of
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.12
Paola Favati, Ornella Menchi
E-ISSN: 2224-2880
77
Volume 21, 2022
the generators of the Toeplitz-like matrix, and finally
in Section 5 the results of the numerical experiments
are shown, confirming the theoretical findings.
Notation: uppercase letters are used for matrix
names, bold lowercase letters are used for vector
names, upper index Tindicates the transpose of a ma-
trix or of a vector, upper index ∗indicates the con-
jugate transpose of a matrix. The special vectors 0
and eiindicate the null vector and the i-th canoni-
cal vector of suitable length. Bold letter idenotes
the imaginary unit. The magnitude of a vector ris
measured by a norm. Specifically krk1=Pi|ri|,
krk2=Pi|ri|2and krk∞=maxi|ri|. Finally, the
symbol denotes the componentwise multiplication
between two vectors of equal length.
2 Circulant and Toeplitz matrices
Circulant matrices and Toeplitz matrices (see [19],
[20], [21]) arise frequently in the numerical treatment
of problems of scientific areas such as physics, signal
and image processing, probability, statistics and many
others. Let us outline some of their properties.
A circulant matrix Mof order kis a square k×k
matrix in which the first row [m1,1, m1,2, . . . , m1,k]T
is given and each subsequent row is rotated one ele-
ment to the right relative to the preceding row. For-
mally, the (i, j)-th element of the i-th row with i=
2, . . . , k is mi,j =mi−1,j−1for j= 2, . . . , k and
mi,1=mi−1,k.
The most important feature of circulant matrices is
that they are diagonalized by a discrete Fourier trans-
form, as shown in the following lemma.
Lemma 1 A circulant matrix Mof size kis diagonal-
ized by the Fourier matrix Fk, whose elements are
fi,j =1
√kω(i−1)(j−1), i, j = 1, . . . , k,
with ω=exp(2πi/k).
Denote by mthe transpose of the first row of M, and
by diag(Fkm)the diagonal matrix whose i-th prin-
cipal element is the i-th element of the vector Fkm.
Then
M=√kFkdiag(Fkm)F∗
k,
and for any vector vit holds
Mv=√kFkFkm F∗
kv.(3)
For a proof of this Lemma, see [19]. 2
The products by Fkand F∗
kcan be efficiently com-
puted by calling FFT with computational cost of order
O(klog k)for k→ ∞.
A Toeplitz matrix Tis a k×kmatrix in which
each descending diagonal from left to right is con-
stant. Two kvectors rand s, with r1=s1, are
given; rTis assumed as first row of Tand sis as-
sumed as first column of T. The (i, j)-th element of
Tis ti,j =t1+j−ifor j≥iand ti,j =t1+i−jfor
i < j. Circulant matrices are special cases of Toeplitz
matrices with s= [r1, rk, . . . , r2]T.
Since the case of triangular Toepliz matrices is of
particular interest, we use the following notation. For
a given vector s,L(s)denotes the lower triangular
Toeplitz matrix whose first column is sand U(r)de-
notes the upper triangular Toeplitz matrix whose first
row is rT, as shown in Figure 1. The matrix-vector
L(s) = "s1
s2s1
.
.
..
.
....
sksk−1··· s1#, U(r) = "r1r2··· rk
r1··· rk−1
....
.
.
r1#
Figure 1: Lower and upper triangular Toeplitz matri-
ces of size k
product of a Toeplitz triangular matrix can be com-
puted by exploiting the properties of circulant matri-
ces. To see how relation (3) is exploited, consider first
the case of a lower triangular Toeplitz matrix of size
n. Let L=L(s), with s= [s1, s2, . . . , sn]T. Given
a vector v, let y=Lvbe the vector to be computed.
The vector vis embedded in a 2n-vector b
vand the
matrix Lis embedded in a 2n×2ncirculant matrix M,
whose first row is the 2n-vector [s1,0T, sn, . . . , s2],
with
b
v=v
0, M =Lb
L
b
L L ,
where b
Lturns out to be an upper triangular Toeplitz
matrix. Since
w=Mb
v=Lv
b
Lv,
the vector yis found in the first half of the vector w.
Then, using (3) the product ycan be computed by the
function lowert given in Algorithm 1.
Algorithm 1: product of a lower triangular
Toeplitz matrix L(s)by a vector v
function lowert (n, s,v)
m= [s1,0T, sn, . . . , s2]T;b
v=v
0;
z=F2nm;p=F∗
2nb
v;
q=zp;w=√2nF2nq;
return w(1 : n);
A similar procedure applies to the upper triangu-
lar Toeplitz matrix U=U(r)whose first row is
rT= [r1, r2, . . . , rn]. Given a vector v, let y=Uv
be the vector to be computed. Proceeding with the
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.12
Paola Favati, Ornella Menchi
E-ISSN: 2224-2880
78
Volume 21, 2022
embedding as before, we see that the circulant matrix
Mhas first row [r1, . . . , rn,0T]. So ycan be com-
puted by the function uppert given in Algorithm 2.
Algorithm 2: product of an upper triangular
Toeplitz matrix U(r)by a vector v
function uppert (n, r,v)
m=r
0;b
v=v
0;
z=F2nm;p=F∗
2nb
v;
q=zp;w=√2nF2nq;
return w(1 : n);
Since the Toeplitz structure is not maintained when
simple operations like multiplication or inversion are
applied, some generalizations have been proposed to
deal with this aspect. Among them, we consider here
a Toeplitz-like structure.
3 Toeplitz-like matrices
The definition of Toeplitz-like structure is based on
the concept of displacement rank [1, 2, 3, 4], which
depends on a particularly chosen displacement oper-
ator and measures how close a matrix is to a Toeplitz
matrix. Given an n×nmatrix A, we consider here
the displacement operator
∇(A) = A−ZAZT,(4)
where Zis the n×ndown-shift matrix, i.e. the binary
matrix with ones only on the subdiagonal and zeros
elsewhere.
The matrix Ais said to be Toeplitz-like if the quan-
tity rdisp(A) = rank ∇(A)(called displacement rank)
is small with respect to n(more formally rdisp(A) =
O(1) for n→ ∞). Let ρ=rdisp(A), then
∇(A) = C DT,(5)
for suitable n×ρmatrices Cand D, called genera-
tors of A. Denoting by ciand di,i= 1, . . . , ρ, the
columns of Cand Drespectively, then
∇(A) =
ρ
X
i=1
cidT
i.(6)
In this sense, we can say that Ais represented through
the generators ci,di. In particular, a Toeplitz matrix
Aof elements ai,j with a1,16= 0, is so represented
∇(A) = c1eT
1+e1dT
2,
with c1=Ae1,d2=ATe1−a11e1,
i. e. c1is the first column of Aand dT
2is the first
row of Awith the first component set to zero. This
shows that the displacement rank of a Toeplitz matrix
is ρ= 2, except in the case of a triangular matrix
where ρ= 1.
The set of Toeplitz-like matrices, unlike the set of
Toeplitz matrices, is closed for the most common op-
erations applied in the numerical algorithms. This
does not mean that the displacement rank is main-
tained during the computation. For example, the
matrix obtained by multiplying two matrices having
rdisp =ρhas the same rdisp, while the displacement
rank of the inverse of a matrix which has rdisp =ρ
may rise up to ρ+ 2. The leading principal submatrix
of a matrix which has rdisp =ρand its Schur comple-
ment have still rdisp =ρ(see [2, 7]).
The generators enable us to express a Toeplitz-like
matrix as the sum of products of lower and upper tri-
angular Toeplitz factors. In fact, from (4) we have
A=∇(A) + ZAZT=∇(A) + Z∇(A)ZT
+Z2A(Z2)T=. . . =
n−1
X
s=0
Zs∇(A)(Zs)T.
Since
n−1
X
s=0 ZsciZsdiT=L(ci)U(di),from (6)
it follows that
A=
ρ
X
i=1
L(ci)U(di),
and, given an n-vector v, we have
Av=
ρ
X
i=1
L(ci)U(di)v.(7)
Using (7) we can compute Avby calling alternatively
the matrix-vector products of upper and lower trian-
gular Toeplitz matrices, as outlined in Algorithm 3. A
saving of the cost can be achieved by skipping the last
FFT call of lowert and exploiting the linearity of F2n
in the final sum.
If the matrix Ais known to be Toeplitz-like, but
it is only given explicitly, any factorization of ∇(A)
can be employed to detect ρand to construct the gen-
erators. For example, we can compute the Gaussian
factorization ∇(A) = LU, where Land Uare lower
and upper triangular matrices, employing a diagonal
pivoting strategy and stopping at the first null pivot.
It follows that while ρ=rank ∇(A)is uniquely
determined, the decomposition (5) of ∇(A), and con-
sequently the representation of Athrough the genera-
tors, is not unique. An important representation is the
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.12
Paola Favati, Ornella Menchi
E-ISSN: 2224-2880
79
Volume 21, 2022
Algorithm 3: product of a Toeplitz-like matrix
by a vector
function prod (n, ρ, C, D, v);
for i= 1 to ρ;
hi=uppert (n, di,v);
gi=lowert (n, ci,hi);
end for;
return
ρ
P
i=1
gi;
orthogonal one [22], obtained by computing the SVD
decomposition ∇(A) = UΣV T, where ρis the num-
ber of strictly positive singular values σiin Σand the
matrices Uand Vhave orthogonal columns.
Denoting by b
Σthe n×ρmatrix having the i-th
principal element equal to √σifor i= 1, . . . , ρ, and
zero elsewere, the n×ρmatrices
Cort =Ub
Σ, Dort =Vb
Σ, (8)
have orthogonal columns and can be assumed as the
generators Cand Din (5).
The following relations hold between the magni-
tudes of Aand ∇(A)
k∇(A)k2≤2kAk2,and kAk2≤nk∇(A)k2.(9)
To measure the magnitude of Awhen it is represented
through the generators, we consider the function
ψ(C, D) =
ρ
X
i=1 kcidT
ik2=
ρ
X
i=1 kcik2kdik2,(10)
which obviously verifies k∇(A)k2≤ψ(C, D).
It is known that the stability of a method solving
a linear system depends on the growth of the matrix
factors computed by the method. If the computations
are performed on Toeplitz-like matrices represented
through the generators, we expect the stability to de-
pend on how large the generators become [23]. In the
general case, no upper bound of ψ(C, D)in terms of
kAk2can be given. However, if the decomposition
(5) is orthogonal, then from (8)
k∇(A)k2=σ1,ci=√σiui,di=√σivi,
where uiis the i-th column of Uand viis the i-th
column of V. Then
kcidT
ik2=σikuik2kvik2=σi,
hence from (9)
ψ(C, D) =
ρ
X
i=1
σi≤ρk∇(A)k2≤2ρkAk2.(11)
4 Stability of the function prod
For the stability analysis we assume that the compu-
tations are carried out in a floating point arithmetic
with unit roundoff . The computed value of a vari-
able (scalar, vector or matrix) vwill be denoted by
evor by “fl(v)”. We assume also that the quanti-
ties which appear in the bounds are not so large to
invalidate a first order error analysis. For simplic-
ity the term “+O(2)”, which appears in the thesis of
the theorems, is omitted in the proofs. Consequently,
any expression of the form xey, where x=O()and
ey−y=O(), is replaced by x y.
The following bounds are used [18]:
•For any vector sit is kL(s)k1=kL(s)k∞=
ksk1, then
kL(s)k2≤pkL(s)k1kL(s)k∞
=ksk1≤√nksk2,(12)
and analogously kU(r)k2≤√nkrk2for any vector
r.
•Given a vector x, a vector whose components
are bounded in modulus by exists such that
x=e
x−e
x+O(2).(13)
•Given two vectors xand y, with e
x=x+δxand
e
y=y+δy, then
fl e
xe
y=xy+θ+O(2),(14)
where
θ=xδy+δxy+xy,
and is a vector whose components are bounded in
modulus by .
•Given ρscalars αiand ρvectors xi,i= 1, . . . , ρ,
then ρvectors χi,i= 1, . . . , ρ, with entries bounded
in modulus by ρ , exist such that
fl ρ
P
i=1
αixi
=
ρ
P
i=1
αixi+xiχi+O(2).
(15)
The following stability result applies to FFT [17]:
•Given a 2n-vector x, let y=F2nxand e
y=
flF2nx, then a 2n×2nmatrix Φexists such that
e
y=y+Φy+O(2),(16)
with kΦk2≤10.7log2(2n). An analogous bound
holds for F∗
2n, with Φreplaced by a matrix Φ∗, which
satisfies the same bound.
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.12
Paola Favati, Ornella Menchi
E-ISSN: 2224-2880
80
Volume 21, 2022
Theorem 2 shows how the computed product of a
triangular Toeplitz matrix by a vector can be regarded
as the exact product of a slightly perturbed matrix by
the vector.
Theorem 2 Given two n-vectors rand v, let U(r)be
the n×nupper triangular Toeplitz matrix whose first
row is rT,
u=U(r)vand e
u=fluppert(n, r,v).
Then a matrix H(r)exists such that
e
u=u+H(r)v+O(2),
where
kH(r)k2≤ γ 0krk2,(17)
with γ0= 42.5√nlog2(2n).
Proof. Applying algorithm uppert we get
e
z=fl F2nm,e
p=fl F∗
2nb
v,
e
q=fl e
ze
p,e
w=fl √2nF2ne
q.
Using (14) and (16) we have
e
z=z+Φz,e
p=p+Φ∗p,
e
q=zp+θ,e
w=√2nF2ne
q+Φw,
where
θ=Φzp+zΦ∗p+zp=Λp,
with
Λ=diag Φz+diag zΦ∗+diag z.
Then
e
w=√2nF2nq+√2nF2nθ+Φw
=w+√2nF2nΛp+Φw.
The vector e
uis found in the first half of e
w. Denoting
by Ethe first half of the identity matrix of order 2n,
we have
e
u=ETe
w,p=F∗
2nEvand e
u=u+H(r)v
with H(r) = √2nETF2nΛF∗
2nE+ETΦM,
where Mis the circulant matrix whose first row is
[rT,0T]. Using (12) and (16) we get
kH(r)k2≤√2nkΛk2+√nkΦk2krk2
≤√2nkΦk2+kΦ∗k2++√nkΦk2krk2
≤42.5√nlog2(2n)krk2.2
An analogous result holds for the matrix-vector
product of a lower triangular Toeplitz computed by
applying algorithm lowert.
Theorem 3 shows how the matrix-vector product
of a Toeplitz-like matrix, computed by the function
prod of Section 3, can be regarded as the exact prod-
uct of a slightly perturbed Toeplitz-like matrix by the
vector.
Theorem 3 Given an n×nToeplitz-like matrix A
with ∇(A) = C DTand an n-vector v, let
u=Avand e
u=flprod(n, ρ, C, D, v).
Then a matrix Θexists such that
e
u=u+Θv+O(2)with kΘk2≤ γ 00 ψ(C, D),
where γ00 =cnlog2(2n),cnot depending on n.
Proof. From (7) we have
u=
ρ
X
i=1
gi,
where for i= 1, . . . , ρ,
gi=L(ci)hiand hi=U(di)v.
The following quantities are effectively computed
e
hi=fluppert (n, di,v),
e
gi=fllowert (n, ci,e
hi,
e
u=flρ
P
i=1 e
gi.
By Theorem 2 we have
e
hi=hi+H(di)v,
where kH(di)k2≤ γ 0kdik2,
e
gi=L(ci)e
hi+H(ci)hi,
where kH(ci)k2≤ γ 0kcik2,then
e
gi=gi+δiv,
where δi=L(ci)H(di) + H(ci)U(di). Using (12)
and (17) we have
kL(ci)k2≤√nkcik2
and
kδik2≤ kL(ci)k2kH(di)k2+kH(ci)k2kU(di)k2
≤2√n γ 0kcik2kdik2.
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.12
Paola Favati, Ornella Menchi
E-ISSN: 2224-2880
81
Volume 21, 2022
Summing for i= 1 , . . . , ρ and applying (15) we get
e
u=flρ
P
i=1 e
gi
=
ρ
P
i=1
gi+
ρ
P
i=1 δiv+giχi=u+Θv,
where the entries of χiare bounded in modulus by ρ
and
Θ=
ρ
X
i=1 δi+diag(χi)L(ci)U(di).
Then from (10)
kΘk2≤2√n γ 0+ρ nρ
P
i=1 kcik2kdik2
≤ γ 00ψ(C, D),where γ00 = 2 √n γ 0+ρ n.
Taking into account the expression of γ0in (17), the
thesis follows. 2
If the decomposition of ∇(A)is orthogonal, then
from (11) it follows
kΘk2≤ γ 00ρk∇(A)k2≤2 γ 00ρkAk2,(18)
suggesting the stability of the algorithm prod.
5 Numerical experiments
The experiments, which have been conducted on an
Intel Core Duo @ 3 GHz, 2GB RAM, using dou-
ble precision arithmetic, have been carried out on
Toeplitz-like matrices of growing size nand different
displacement rank ρ.
Two sets of numerical experiments are performed,
in order to validate the upper bound given in Theorem
3 by investigating the behaviour of the error produced
in the computation of prod (n, ρ, C, D, v).
(i) The matrices for the first set of experiments have
been generated for different values of the displace-
ment rank and growing values of nin the range
[23,29]. For each size nand fixed values of ρ, ten
triples {C, D, v}, with entries uniformly distributed
in [−10,10] and v6=0, are randomly generated. In
Figure 2 the arithmetic mean µnof the errors ke
u−
uk2/kvk2is plotted versus nfor the case ρ= 5 (no
significant differences occurring for other values of
ρ), together with the upper bound τn= γ 00ψ(C, D)
of Theorem 3, with γ00 = 85 nlog2(2n) + 5 n. As
expected, µnis largely overestimated by τn.
(ii) For the second set of experiments we fix n= 29
and ρ= 5 and generate matrices Cand Das in the
previous case, except for the fact that different pairs
of generators corresponding to the same matrix Aare
100
200
300
400
500
-11
-10
-9
-8
-7
-6
Figure 2: Log plot of µn(dashed line) and τn(solid
line) as functions of n.
β ψ(Cβ, Dβ)eβ
101260 5.5 10−11
102391 6.5 10−11
1031.7 1036.1 10−10
1041.5 1046.7 10−8
1051.5 1057.4 10−6
1061.5 1065.6 10−5
1071.5 1079.1 10−2
1081.5 1085.4 100
ort 247 1.2 10−10
Table 1: Values of ψ(Cβ, Dβ)and of eβvarying β.
Last row shows the corresponding values for the or-
thogonal representation.
obtained by allowing the columns of Cand Dto de-
pend on a parameter β. In this way, very different val-
ues of the function ψ(Cβ, Dβ)occur, which increase
with β. Table 1 shows that also the relative errors
eβ=ke
u−uk2/kvk2increase with β. However by
using the orthogonal representation Cort,Dort of A,
the function ψ(Cort, Dort)can be bounded by kAk2
(in this example kAk2= 355) and consequently the
relative error is reduced as suggested by (18) (see last
row of Table 1).
6 Conclusions
The numerical stability of the matrix-vector product
for Toeplitz-like matrices, performed via FFT, has
been analyzed. The analysis has pointed out that the
error greatly depends on the magnitude of the gener-
ators of the matrix. The numerical experimentation
confirms this result, suggesting that the magnitude of
the generators should be monitored, and the genera-
tors should be replaced by orthogonal ones when they
become too large with respect to the magnitude of
the associated matrix. For example, when ψ(C, D)
becomes larger than 2ρkAk2. The present study is
part of a research which analyzes some superfast algo-
rithms, like the one proposed in [7], for the solution of
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.12
Paola Favati, Ornella Menchi
E-ISSN: 2224-2880
82
Volume 21, 2022
Toeplitz-like systems from the stability point of view.
From a theoretical error analysis, it turns out that in-
stability and computational complexity balance, since
in general larger errors tend to be produced by faster
algorithms. For the near future we are planning to
continue our study of this interesting field, by detect-
ing the parameters which rule the stability behavior
of iterative superfast algorithms and focusing on the
conditioning properties connected with the magnitude
of the generators of the matrices involved.
References:
[1] T. Kailath, A. Viera and M. Morf, “Inverses of
Toeplitz operators, innovations and orthogonal
polynomials”, SIAM Rev., 20, pp. 106-119, 1978.
[2] T. Kailath, S. Y. Kung and M. Morf, “Displace-
ment ranks of matrices and linear equations”, J.
Math. Anal. Appl., 68, pp. 395-407, 1979.
[3] G. Heinig and K. Rost, Algebraic methods for
Toeplitz-like matrices and operators, Akademie-
Verlag, Berlin, 1984.
[4] T. Kailath and A. H. Sayed, “Displacement struc-
ture: theory and applications”, SIAM Rev., 37, pp.
297-386, 1995.
[5] D. A. Bini, “Matrix structures and applica-
tions”, Centre international de rencontres mathe-
matiques, U.M.S. 822 C.N.R.S./S.M.F., 4, pp. 1-
45, 2014.
[6] A. Aricó and G. Rodriguez, “ A fast solver for lin-
ear systems with displacement structure”, Numer.
Algorithms, 55, pp. 529-556, 2010.
[7] P. Favati, G. Lotti and O. Menchi, “A Divide
and Conquer Algorithm for the Superfast solution
of Toeplitz-like Systems”, SIAM. J. Matrix Anal.
Appl., 33, pp. 1039-1056, 2012.
[8] Y. Xi, J. Xia, S.Cauley andV.Balakrishnan, “Su-
perfast and stable structured solvers for Toeplitz
least squares via randomized sampling”, SIAM. J.
Matrix Anal. Appl., 35, pp. 44-72, 2014.
[9] R. Ke, M. K. Ng and H. W. Sun, “A fast direct
method for block triangular Toeplitz-like with tri-
diagonal block systems from time-fractional par-
tial differential equations”, Journal of Computa-
tional Physics, 303, pp. 203-211, 2015.
[10] X. Lin, M. K. Ng and H. W. Sun, “A Split-
ting Preconditioner for Toeplitz-Like Linear Sys-
tems Arising from Fractional Diffusion Equa-
tions”, SIAM. J. Matrix Anal. Appl., 38, pp. 1580-
1614, 2017.
[11] N. Akhoundi, “Toeplitz-like preconditioner for
linear systems from spatial fractional diffusion
equations”, Iranian Journal of Numerical Anal-
ysis and Optimization, 11, pp. 95-106, 2021.
[12] D. Bini and B. Meini, “On Cyclic Reduction Ap-
plied to a Class of Toeplitz-Like Matrices Arising
in Queueing Problems”, in: W. J. Stewart (eds)
Computations with Markov Chains. Springer,
Boston, MA , 1995.
[13] A. Böttcher, C. Garoni and S. Serra-Capizzano,
“Exploration of Toeplitz-like matrices with un-
bounded symbols is not a purely academic jour-
ney”, Sb. Math., 208, pp.1602-1627, 2017.
[14] A. Bostan, C. P. Jeannerod, C. Mouilleron
and E. Schost, “On matrices with displacement
structure: generalized operators and faster algo-
rithms”, SIAM. J. Matrix Anal. Appl., 38, pp. 733-
775, 2017.
[15] G. J. Groenewald, S. ter Horst, J. Jaftha and A.
C. M. Ran, “A Toeplitz-Like Operator with Ra-
tional Matrix Symbol Having Poles on the Unit
Circle: Fredholm Properties”, Complex Analysis
and Operator Theory 15, pp. 1-29, 2021.
[16] J. W. Cooley and O. W. Tukey, “An Algorithm
for the Machine Calculation of Complex Fourier
Series”, Math. Comput., 19, pp. 297-301, 1965.
[17] M. Arioli, H. Munthe-Kaas and L. Valdettaro,
“Componentwise error analysis for FFT’s with
applications to fast Helmholtz solvers”, Numer.
Algorithms, 12, pp. 65-88, 1996.
[18] N. J. Higham, Accuracy and Stability of Numer-
ical Algorithms, SIAM, Philadelphia, PA, 1996.
[19] P. J. Davis, Circulant Matrices, Wiley, New
York, 1979.
[20] G. H. Golub and C. F. Van Loan, “Circulant
Systems”, par. 4.7.7 in Matrix Computations (3rd
ed.), Johns Hopkins, 1996.
[21] R. M. Gray, “Toeplitz and Circulant Matrices: A
Review”, Foundations and Trends in Communi-
cations and Information Theory, 2, pp. 155-239,
2006.
[22] V. Y. Pan, Y. Rami and X. Wang, “Structured ma-
trices and Newton’s iteration: unified approach”,
Linear Algebra and its Applications, 343-344, pp.
233-265, 2002.
[23] P. Favati, G. Lotti and O. Menchi, “Stability
of the Levinson algorithm for Toeplitz-like sys-
tems”, SIAM Journal on Matrix Analysis and Ap-
plications, 31, pp. 2531-2552, 2010.
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.12
Paola Favati, Ornella Menchi
E-ISSN: 2224-2880
83
Volume 21, 2022
Contribution of individual authors to
the creation of a scientific article
(ghostwriting policy)
All the authors have contributed substantially and in
equal measure to all the phases of the work reported.
Creative Commons Attribution
License 4.0 (Attribution 4.0
International , CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/li-
censes/by/4.0/deed.en_US
WSEAS TRANSACTIONS on MATHEMATICS
DOI: 10.37394/23206.2022.21.12
Paola Favati, Ornella Menchi
E-ISSN: 2224-2880
84
Volume 21, 2022
