V-cycle optimal convergence for DCT-III matrices
7 V-cycle optimal convergence for DCT-III ma-
trices
C. Tablino Possio
Dedicated to Georg Heinig
Abstract. The paper analyzes a two-grid and a multigrid method for matrices
belonging to the DCT-III algebra and generated by a polynomial symbol. The
aim is to prove that the convergence rate of the considered multigrid method
(V-cycle) is constant independent of the size of the given matrix. Numerical
examples from differential and integral equations are considered to illustrate
the claimed convergence properties.
Mathematics Subject Classification (2000). Primary 65F10, 65F15, 15A12.
Keywords. DCT-III algebra, two-grid and multigrid iterations, multi-iterative
methods.
1. Introduction
In the last two decades, an intensive work has concerned the numerical solution
of structured linear systems of large dimensions [6, 14, 16]. Many problems have
been solved mainly by the use of (preconditioned) iterative solvers. However, in
the multilevel setting, it has been proved that the most popular matrix algebra
preconditioners cannot work in general (see [23, 26, 20] and references therein).
On the other hand, the multilevel structures often are the most interesting in
practical applications. Therefore, quite recently, more attention has been focused
(see [1, 2, 7, 5, 27, 9, 12, 10, 13, 22, 25, 19]) on the multigrid solution of multilevel
structured (Toeplitz, circulants, Hartley, sine (τ class) and cosine algebras) linear
systems in which the coefficient matrix is banded in a multilevel sense and positive
definite. The reason is due to the fact that these techniques are very efficient, the
total cost for reaching the solution within a preassigned accuracy being linear as
the dimensions of the involved linear systems.
The work of the author was partially supported by MIUR, grant number 2006017542.
http://arxiv.org/abs/0704.1980v1
2 C. Tablino Possio
In this paper we deal with the case of matrices generated by a polynomial symbol
and belonging to the DCT-III algebra. This kind of matrices appears in the solution
of differential equations and integral equations, see for instance [4, 18, 24]. In
particular, they directly arise in certain image restoration problems or can be used
as preconditioners for more complicated problems in the same field of application
[17, 18].
In [7] a Two-Grid (TGM)/Multi-Grid (MGM) Method has been proposed
and the theoretical analysis of the TGM has been performed in terms of the alge-
braic multigrid theory developed by Ruge and Stüben [21].
Here, the aim is to provide general conditions under which the proposed MGM
results to be optimally convergent with a convergence rate independent of the di-
mension and to perform the corresponding theoretical analysis.
More precisely, for MGM we mean the simplest (and less expensive) version of the
large family of multigrid methods, i.e. the V-cycle procedure. For a brief descrip-
tion of the TGM and of the MGM (standard V-cycle) we refer to §2. An extensive
treatment can be found in [11], and especially in [28].
In all the considered cases the MGM results to be optimal in the sense of Defi-
nition 1.1, i.e. the problem of solving a linear system with coefficient matrix Am
is asymptotically of the same cost as the direct problem of multiplying Am by a
vector.
Definition 1.1. [3] Let {Amxm = bm} be a given sequence of linear systems of
increasing dimensions. An iterative method is optimal if
1. the arithmetic cost of each iteration is at most proportional to the complexity
of a matrix vector product with matrix Am,
2. the number of iterations for reaching the solution within a fixed accuracy can
be bounded from above by a constant independent of m.
In fact, the total cost of the proposed MGM will be of O(m) operations since for
any coarse level s we can find a projection operator P ss+1 such that
• the matrix vector product involving P ss+1 costs O(ms) operations wherems =
m/2s;
• the coarse grid matrix Ams+1 = P ss+1Ams(P ss+1)T is also a matrix in the
DCT III algebra generated by a polynomial symbol and can be formed within
O(ms) operations;
• the convergence rate of the MGM is independent of m.
The paper is organized as follows. In §2 we briefly report the main tools re-
garding to the convergence theory of algebraic multigrid methods [21]. In §3 we
consider the TGM for matrices belonging to DCT-III algebra with reference to
some optimal convergence properties, while §4 is devoted to the convergence anal-
ysis of its natural extension as V-cycle. In §5 numerical evidences of the claimed
results are discussed and §6 deals with complexity issues and conclusions.
V-cycle optimal convergence for DCT-III matrices 3
2. Two-grid and Multi-grid methods
In this section we briefly report the main results pertaining to the convergence
theory of algebraic multigrid methods.
Let us consider the generic linear system Amxm = bm, where Am ∈ Cm×m is a
Hermitian positive definite matrix and xm, bm ∈ Cm. Let m0 = m > m1 > . . . >
ms > . . . > msmin and let P
s+1 ∈ Cms+1×ms be a given full-rank matrix for any
s. Lastly, let us denote by Vs a class of iterative methods for linear systems of
dimension ms.
According to [11], the algebraic Two-Grid Method (TGM) is an iterative method
whose generic step is defined as follow.
xouts = T GM(s, xins , bs)
xpres = V
s,pre(x
s ) Pre-smoothing iterations
rs = Asx
s − bs
rs+1 = P
s+1rs
As+1 = P
s+1As(P
Solve As+1ys+1 = rs+1
x̂s = x
s − (P ss+1)Hys+1
Exact Coarse Grid Correction
xouts = V
νpost
s,post(x̂s) Post-smoothing iterations
where the dimension ms is denoted in short by the subscript s.
In the first and last steps a pre-smoothing iteration and a post-smoothing iteration
are respectively applied νpre times and νpost times, according to the chosen iterative
method in the class Vs. Moreover, the intermediate steps define the so called exact
coarse grid correction operator, that depends on the considered projector operator
P ss+1. The global iteration matrix of the TGM is then given by
TGMs = V
νpost
s,postCGCsV
s,pre, (2.1)
CGCs = Is − (P ss+1)HA−1s+1P
s+1As As+1 = P
s+1As(P
H , (2.2)
where Vs,pre and Vs,post respectively denote the pre-smoothing and post-smoothing
iteration matrices.
4 C. Tablino Possio
By means of a recursive procedure, the TGM gives rise to a Multi-Grid
Method (MGM): the standard V-cycle is defined as follows.
xouts = MGM(s, xins , bs)
if s ≤ smin then
Solve Asx
s = bs Exact solution
xpres = V
s,pre(x
s ) Pre-smoothing iterations
rs = Asx
s − bs
rs+1 = P
s+1rs
ys+1 = MGM(s+ 1,0s+1, rs+1)
x̂s = x
s − (P ss+1)Hys+1
Coarse Grid Correction
xouts = V
νpost
s,post(x̂s) Post-smoothing iterations
Notice that in MGM the matrices As+1 = P
s+1As(P
H are more profitably
formed in the so called setup phase in order to reduce the computational costs.
The global iteration matrix of the MGM can be recursively defined as
MGMsmin = O ∈ Csmin×smin ,
MGMs = V
νpost
s,post
Is − (P ss+1)H (Is+1 −MGMs+1)A−1s+1P ss+1As
s,pre,
s = smin − 1, . . . , 0.
Some general conditions that ensure the convergence of an algebraic TGM
and MGM are due to Ruge and Stüben [21].
Hereafter, by ‖·‖2 we denote the Euclidean norm on Cm and the associated induced
matrix norm over Cm×m. If X is positive definite, ‖ · ‖X = ‖X1/2 · ‖2 denotes the
Euclidean norm weighted by X on Cm and the associated induced matrix norm.
Finally, if X and Y are Hermitian matrices, then the notation X ≤ Y means that
Y −X is nonnegative definite.
Theorem 2.1 (TGM convergence [21]). Let m0, m1 be integers such that m0 >
m1 > 0, let A ∈ Cm0×m0 be a positive definite matrix. Let V0 be a class of iterative
methods for linear systems of dimension m0 and let P
1 ∈ Cm1×m0 be a given full-
rank matrix. Suppose that there exist αpre > 0 and αpost > 0 independent of m0
V-cycle optimal convergence for DCT-III matrices 5
such that
‖V0,pre x‖2A ≤ ‖x‖2A − αpre‖V0,pre x‖2AD−1A for any x ∈ C
m0 (2.3a)
‖V0,post x‖2A ≤ ‖x‖2A − αpost ‖x‖2AD−1A for any x ∈ C
m0 (2.3b)
(where D denotes the main diagonal of A) and that there exists γ > 0 independent
of m0 such that
y∈Cm1
‖x− (P 01 )Hy‖2D ≤ γ‖x‖2A for any x ∈ Cm0 . (2.4)
Then, γ ≥ αpost and
‖TGM0‖A ≤
1− αpost/γ
1 + αpre/γ
. (2.5)
It is worth stressing that in Theorem 2.1 the matrix D ∈ Cm0×m0 can be substi-
tuted by any Hermitian positive definite matrix X : clearly the choice X = I can
give rise to valuable simplifications [1].
At first sight, the MGM convergence requirements are more severe since the
smoothing and CGC iteration matrices are linked in the same inequalities as stated
below.
Theorem 2.2 (MGM convergence [21]). Let m0 = m > m1 > m2 > . . . > ms >
. . . > msmin and let A ∈ Cm×m be a positive definite matrix. Let P ss+1 ∈ Cms+1×ms
be full-rank matrices for any level s. Suppose that there exist δpre > 0 and δpost > 0
such that
‖V νpres,prex‖
≤ ‖x‖2As − δpre ‖CGCsV
s,prex‖
for any x ∈ Cms (2.6a)
‖V νposts,posrx‖
≤ ‖x‖2As − δpost ‖CGCsx‖
for any x ∈ Cms (2.6b)
both for each s = 0, . . . , smin − 1, then δpost ≤ 1 and
‖MGM0‖A 6
1− δpost
1 + δpre
< 1. (2.7)
By virtue of Theorem 2.2, the sequence {x(k)m }k∈N will converge to the solution of
the linear system Amxm = bm and within a constant error reduction not depending
on m and smin if at least one between δpre and δpost is independent of m and smin.
Nevertheless, as also suggested in [21], the inequalities (2.6a) and (2.6b) can
be respectively splitted as
‖V νpres,prex‖
≤ ‖x‖2As − α ‖V
s,prex‖AsD−1s As
‖CGCsx‖2As ≤ γ ‖x‖
δpre = α/γ
(2.8)
and 
‖V νposts,postx‖
≤ ‖x‖2As − β ‖x‖
‖CGCs x‖2As ≤ γ ‖x‖
δpost = β/γ
(2.9)
6 C. Tablino Possio
where Ds is the diagonal part of As (again, the AD
−1A-norm is not compulsory [1]
and the A2-norm will be considered in the following) and where, more importantly,
the coefficients α, β and γ can differ in each recursion level s since the step from
(2.8) to (2.6a) and from (2.9) to (2.6b) are purely algebraic and do not affect the
proof of Theorem 2.2.
Therefore, in order to prove the V-cycle optimal convergence, it is possible to
consider the inequalities
‖V νpres,prex‖
≤ ‖x‖2As − αs ‖V
s,prex‖
for any x ∈ Cms (2.10a)
‖V νposts,postx‖
≤ ‖x‖2As − βs ‖x‖
for any x ∈ Cms (2.10b)
‖CGCsx‖2As ≤ γs ‖x‖
for any x ∈ Cms . (2.10c)
where it is required that αs, βs, γs ≥ 0 for each s = 0, . . . , smin − 1 and
δpre = min
0≤s<smin
, δpost = min
0≤s<smin
. (2.11)
We refer to (2.10a) as the pre-smoothing property, (2.10b) as the post-smoothing
property and (2.10c) as the approximation property (see [21]).
An evident benefit in considering the inequalities (2.10a)-(2.10c) relies on to the
fact that the analysis of the smoothing iterations is distinguished from the more
difficult analysis of the projector operator.
Moreover, the MGM smoothing properties (2.10a) and (2.10b) are nothing more
than the TGM smoothing properties (2.3a) and (2.3b) with D substituted by I,
in accordance with the previous reasoning (see [1]).
3. Two-grid and Multi-grid methods for DCT III matrices
Let Cm = {Cm ∈ Rm×m|Cm = QmDmQTm} the unilevel DCT-III cosine matrix
algebra, i.e. the algebra of matrices that are simultaneously diagonalized by the
orthogonal transform
2− δj,1
(i − 1)(j − 1/2)π
i,j=1
(3.1)
with δi,j denoting the Kronecker symbol.
Let f be a real-valued even trigonometric polynomial of degree k and period 2π.
Then, the DCT III matrix of order m generated by f is defined as
Cm(f) = QmDm(f)Q
m, Dm(f) = diag1≤j≤m f
(j − 1)π
Clearly, Cm(f) is a symmetric band matrix of bandwidth 2k+1. In the following,
we denote in short with Cs = Cms(gs) the DCT III matrix of size ms generated
by the function gs.
An algebraic TGM/MGM method for (multilevel) DCT III matrices generated by
a real-valued even trigonometric polynomial has been proposed in [7]. Here, we
V-cycle optimal convergence for DCT-III matrices 7
briefly report the relevant results with respect to TGM convergence analysis, the
aim being to prove in §4 the V-cycle optimal convergence under suitable conditions.
Indeed, the projector operator P ss+1 is chosen as
P ss+1 = T
s+1Cs(ps)
where T ss+1 ∈ Rms+1×ms , ms+1 = ms/2, is the cutting operator defined as
T ss+1
2 for j ∈ {2i− 1, 2i}, i = 1, . . . ,ms+1,
0 otherwise.
(3.2)
and Cs(ps) is the DCT-III cosine matrix of size ms generated by a suitable even
trigonometric polynomial ps.
Here, the scaling by a factor 1/
2 is introduced in order to normalize the matrix
T ss+1 with respect to the Euclidean norm. From the point of view of an algebraic
multigrid this is a natural choice, while in a geometric multigrid it is more natural
to consider just a scaling by 1/2 in the projector, to obtain an average value.
The cutting operator plays a leading role in preserving both the structural and
spectral properties of the projected matrix Cs+1: in fact, it ensures a spectral link
between the space of the frequencies of size ms and the corresponding space of
frequencies of size ms+1, according to the following Lemma.
Lemma 3.1. [7] Let Qs ∈ Rms×ms and T ss+1 ∈ Rms+1×ms be given as in (3.1) and
(3.2) respectively. Then
T ss+1Qs = Qs+1[Φs+1,Θs+1Πs+1], (3.3)
where
Φs+1 = diagj=1,...,ms+1
(j − 1)π
, (3.4a)
Θs+1 = diagj=1,...,ms+1
− cos
, (3.4b)
and Πs+1 ∈ Rms+1×ms+1 is the permutation matrix
(1, 2, . . . ,ms+1) 7→ (1,ms+1,ms+1 − 2, . . . , 2).
As a consequence, let As = Cs(fs) be the DCT-III matrix generated by fs, then
As+1 = P
s+1As(P
T = Cs+1(fs+1)
where
fs+1(x) = cos
(3.5)
+ cos2
π − x/2
π − x
π − x
, x ∈ [0, π].
On the other side, the convergence of proposed TGM at size ms is ensured
by choosing the polynomial as follows.
8 C. Tablino Possio
Definition 3.2. Let x0 ∈ [0, π) a zero of the generating function fs. The polynomial
ps is chosen so that
p2s(π − x)
fs(x)
< +∞, (3.6a)
p2s(x) + p
s(π − x) > 0. (3.6b)
In the special case x0 = π, the requirement (3.6a) is replaced by
x→x0=π
p2s(π − x)
fs(x)
< +∞. (3.7a)
If fs has more than one zero in [0, π], then ps will be the product of the polynomials
satisfying the condition (3.6a) (or (3.7a)) for every single zero and globally the
condition (3.6b).
It is evident from the quoted definition that the polynomial ps must have
zeros of proper order in any mirror point x̂0 = π − x0, where x0 is a zeros of fs.
It is worth stressing that conditions (3.6a) and (3.6b) are in perfect agreement with
the case of other structures such as τ , symmetric Toeplitz and circulant matrices
(see e.g. [22, 25]), while the condition (3.7a) is proper of the DCT III algebra and
it corresponds to a worsening of the convergence requirements.
Moreover, as just suggested in [7], in the case x0 = 0 the condition (3.6a) can also
be weakened as
x→x0=0
p2s(π − x)
fs(x)
< +∞. (3.8a)
We note that if fs is a trigonometric polynomial of degree k, then fs can have a
zero of order at most 2k. If none of the root of fs are at π, then by (3.6a) the
degree of ps has to be less than or equal to ⌈k/2⌉. If π is one of the roots of fs,
then the degree of ps is less than or equal to ⌈(k + 1)/2⌉.
Notice also that from (3.5), it is easy to obtain the Fourier coefficients of fs+1 and
hence the nonzero entries of As+1 = Cs+1(fs+1). In addition, we can obtain the
roots of fs+1 and their orders by knowing the roots of fs and their orders.
Lemma 3.3. [7] If 0 ≤ x0 ≤ π/2 is a zero of fs, then by (3.6a), ps(π − x0) = 0
and hence by (3.5), fs+1(2x
0) = 0, i.e. y0 = 2x0 is a zero of fs+1. Furthermore,
because ps(π − x0) = 0, by (3.6b), ps(x0) > 0 and hence the orders of x0 and y0
are the same. Similarly, π/2 ≤ x0 < π, then y0 = 2(π − x0) is a root of fs+1 with
the same order as x0. Finally, if x0 = π, then y0 = 0 with order equals to the
order of x0 plus two.
In [7] the Richardson method has be considered as the most natural choice
for the smoothing iteration, since the corresponding iteration matrix Vm := Im −
ωAm ∈ Cm×m belongs to the DCT-III algebra, too. Further remarks about such
a type of smoothing iterations and the tuning of the parameter ω are reported in
[25, 2].
V-cycle optimal convergence for DCT-III matrices 9
Theorem 3.4. [7] Let Am0 = Cm0(f0) with f0 being a nonnegative trigonometric
polynomial and let Vm0 = Im0 − ωAm0 with ω = 2/‖f0‖∞ and ω = 1/‖f0‖∞,
respectively for the pre-smoothing and the post-smoothing iteration. Then, under
the quoted assumptions and definitions the inequalities (2.3a), (2.3b), and (2.4)
hold true and the proposed TGM converges linearly.
Here, it could be interesting to come back to some key steps in the proof of
the quoted Theorem 3.4 in order to highlight the structure with respect to any
point and its mirror point according to the considered notations.
By referring to a proof technique developed in [22], the claimed thesis is obtained
by proving that the right-hand sides in the inequalities
γ ≥ 1
ds(x)
π − x
p2s(π − x)
fs(x)
, (3.9a)
γ ≥ 1
ds(x)
π − x
p2s(π − x)
fs(x)
+ cos2
) p2s(x)
fs(π − x)
, (3.9b)
ds(x) = cos
p2s(x) + cos
π − x
p2s(π − x) (3.9c)
are uniformly bounded on the whole domain so that γ is an universal constant.
It is evident that (3.9a) is implied by (3.9b). Moreover, both the two terms in
(3.9b) and in ds(x) can be exchanged each other, up to the change of variable
y = π − x.
Therefore, if x0 6= π it is evident that Definition 3.2 ensures the required uniform
boundedness since the condition p2s(x) + p
s(π − x) > 0 implies ds(x) > 0.
In the case x0 = π, the inequality (3.9b) can be rewritten as
γ ≥ 1
p2s(x)
p2s(π − x)
p2s(π − x)
fs(x)
p2s(x)
fs(π − x)
(3.10)
so motivating the special case in Definition 3.2.
4. V-cycle optimal convergence
In this section we propose a suitable modification of Definition 3.2 with respect to
the choice of the polynomial involved into the projector, that allows us to prove
the V-cycle optimal convergence according to the verification of the inequalities
(2.10a)-(2.10c) and the requirement (2.11).
It is worth stressing that the MGM smoothing properties do not require a true
verification, since (2.10a) and (2.10b) are exactly the TGM smoothing properties
(2.3a) and (2.3b) (with D = I).
Proposition 4.1. Let As = Cms(fs) for any s = 0, . . . , smin, with fs ≥ 0, and
let ωs be such that 0 < ωs ≤ 2/‖fs‖∞. If we choose αs and βs such that αs ≤
10 C. Tablino Possio
ωsmin
2, (2− ωs‖fs‖∞)/(1− ωs‖fs‖∞)2
and βs 6 ωs(2 − ωs‖fs‖∞) then for
any x ∈ Cm the inequalities
‖Vs,pre x‖2As ≤ ‖x‖
− αs ‖Vs,pre x‖2As (4.1)
‖Vs,post x‖2As ≤ ‖x‖
− βs ‖x‖2As (4.2)
hold true.
Notice, for instance, that the best bound to βs is given by 1/‖fs‖∞ and it is
obtained by taking ωs = 1/‖fs‖∞ [25, 2].
Concerning the analysis of the approximation condition (2.10c) we consider
here the case of a generating function f0 with a single zero at x
0. In such a case,
the choice of the polynomial in the projector is more severe with respect to the
case of TGM.
Definition 4.2. Let x0 ∈ [0, π) a zero of the generating function fs. The polynomial
ps is chosen in such a way that
ps(π − x)
fs(x)
< +∞, (4.3a)
p2s(x) + p
s(π − x) > 0. (4.3b)
In the special case x0 = π, the requirement (4.3a) is replaced by
x→x0=π
ps(π − x)
fs(x)
< +∞. (4.4a)
Notice also that in the special case x0 = 0 the requirement (4.3a) can be weakened
x→x0=0
ps(π − x)
fs(x)
< +∞. (4.5a)
Proposition 4.3. Let As = Cms(fs) for any s = 0, . . . , smin, with fs ≥ 0. Let
P ss+1 = T
s+1Cs(ps), where ps(x) is fulfilling (4.3a) (or (4.4a)) and (4.3b). Then,
for any s = 0, . . . , smin − 1, there exists γs > 0 independent of ms such that
‖CGCsx‖2As ≤ γs ‖x‖
for any x ∈ Cms , (4.6)
where CGCs is defined as in (2.2).
Proof. Since
CGCs = Is − (P ss+1)T (P ss+1As(P ss+1)T )−1P ss+1As
is an unitary projector, it holds that CGCTs As CGCs = As CGCs. Therefore, the
target inequality (4.6) can be simplified and symmetrized, giving rise to the matrix
inequality
C̃GCs = Is −A1/2s (P ss+1)T (P ss+1As(P ss+1)T )−1P ss+1A1/2s ≤ γsAs. (4.7)
V-cycle optimal convergence for DCT-III matrices 11
Hence, by invoking Lemma 3.1, QTs C̃GCsQs can be permuted into a 2 × 2 block
diagonal matrix whose jth block, j = 1, . . . ,ms+1, is given by the rank-1 matrix
(see [8] for the analogous τ case)
c2j + s
c2j cjsj
cjsj s
where
cj = cos
p2f(x
j ) sj = − cos
π − x[ms]j
p2f(π − x[ms]j ).
As in the proof of the TGM convergence, due to the continuity of fs and ps, (4.7)
is proven if the right-hand sides in the inequalities
d̃s(x)
π − x
p2sfs(π − x)
fs(x)
(4.8a)
d̃s(x)
π − x
p2sfs(π − x)
fs(x)
+ cos2
) p2sfs(x)
fs(π − x)
(4.8b)
d̃s(x) = cos
p2sfs(x) + cos
π − x
p2sf(π − x) (4.8c)
are uniformly bounded on the whole domain so that γs are universal constants.
Once again, it is evident that (4.8a) is implied by (4.8b). Moreover, both the terms
in (4.8b) and in d̃s(x) can be exchanged each other, up to the change of variable
y = π − x.
Therefore, if x0 6= π, (4.8b) can be rewritten as
d̂s(x)
π − x
p2s(π − x)
f2s (x)
+ cos2
) p2s(x)
f2s (π − x)
(4.9)
where
d̂s(x) = cos
) p2s(x)
fs(π − x)
+ cos2
π − x
p2s(π − x)
fs(x)
so that Definition 4.2 ensures the required uniform boundedness.
In the case x0 = π, the inequality (4.8b) can be rewritten as
p2s(x)
fs(π − x)
p2s(π − x)
fs(x)
p2s(π − x)
f2s (x)
p2s(x)
f2s (π − x)
(4.10)
so motivating the special case in Definition 4.2. �
12 C. Tablino Possio
Remark 4.4. Notice that in the case of pre-smoothing iterations and under the
assumption Vs,pre nonsingular, the approximation condition
‖CGCsV νpres,prex‖
≤ γs ‖V νpres,prex‖
for any x ∈ Cms , (4.11)
is equivalent to the condition, in matrix form, C̃GCs ≤ γsAs obtained in Propo-
sition 4.3.
In Propositions 4.1 and 4.3 we have obtained that for every s (independent
of m = m0) the constants αs, βs, and γs are absolute values not depending on
m = m0, but only depending on the functions fs and ps. Nevertheless, in order to
prove the MGM optimal convergence according to Theorem 2.2, we should verify
at least one between the following inf–min conditions [1]:
δpre = inf
0≤s≤log2(m0)
> 0, δpost = inf
0≤s≤log2(m0)
> 0. (4.12)
First, we consider the inf-min requirement (4.12) by analyzing the case of a gen-
erating function f̃0 with a single zero at x
0 = 0.
It is worth stressing that in such a case the DCT-III matrix Ãm0 = Cm0(f̃0) is sin-
gular since 0 belongs to the set of grid points x
j = (j − 1)π/m0, j = 1, . . . ,m0.
Thus, the matrix Ãm0 is replaced by
Am0 = Cm0(f0) = Cm0(f̃0) + f̃0
with e = [1, . . . , 1]T ∈ Rm0 and where the rank-1 additional term is known as
Strang correction [29]. Equivalently, f̃0 ≥ 0 is replaced by the generating function
f0 = f̃0 + f̃0
1 +2πZ
> 0, (4.13)
where χX is the characteristic function of the set X and w
1 = x
0 = 0.
In Lemma 4.5 is reported the law to which the generating functions are subjected at
the coarser levels. With respect to this target, it is useful to consider the following
factorization result: let f ≥ 0 be a trigonometric polynomial with a single zero at
x0 of order 2q. Then, there exists a positive trigonometric polynomial ψ such that
f(x) = [1− cos(x− x0)]q ψ(x). (4.14)
Notice also that, according to Lemma 3.3, the location of the zero is never shifted
at the subsequent levels.
Lemma 4.5. Let f0(x) = f̃0(x) + c0χ2πZ(x), with f̃0(x) = [1− cos(x)]qψ0(x), q
being a positive integer and ψ0 being a positive trigonometric polynomial and with
c0 = f̃0
. Let ps(x) = [1+cos(x)]
q for any s = 0, . . . , smin− 1. Then, under
the same assumptions of Lemma 3.1, each generating function fs is given by
fs(x) = f̃s(x) + csχ2πZ(x), f̃s(x) = [1− cos(x)]qψs(x).
V-cycle optimal convergence for DCT-III matrices 13
The sequences {ψs} and {cs} are defined as
ψs+1 = Φq,ps(ψs), cs+1 = csp
s(0), s = 0, . . . , smin − 1,
where Φq,p is an operator such that
[Φq,p(ψ)] (x) =
(ϕpψ)
+ (ϕpψ)
π − x
, (4.15)
with ϕ(x) = 1 + cos(x). Moreover, each f̃s is a trigonometric polynomial that
vanishes only at 2πZ with the same order 2q as f̃0.
Proof. The claim is a direct consequence of Lemma 3.1. Moreover, since the func-
tion ψ0 is positive by assumption, the same holds true for each function ψs. �
Hereafter, we make use of the following notations: for a given function f , we
will write Mf = supx |f |, mf = infx |f | and µ∞(f) =Mf/mf .
Now, if x ∈ (0, 2π) we can give an upper bound for the left-hand side R(x) in
(4.9), since it holds that
R(x) =
p2s(x)
f2s (π − x)
p2s(π − x)
f2s (x)
p2s(x)
fs(π − x)
p2s(π − x)
fs(x)
ψ2s (π − x)
ψ2s(x)
ps(x)
ψs(π − x)
ps(π − x)
ψs(x)
≤ Mψs
ps(x) + cos2
ps(π − x)
≤ Mψs
we can consider γs = Mψs/m
. In the case x = 0, since ps(0) = 0, it holds
R(0) = 1/fs(π), so that we have also to require 1/fs(π) ≤ γs. However, since
1/fs(π) ≤Mψs/m2ψs , we take γ
s =Mψs/m
as the best value.
In (2.9), by choosing ω∗s = ‖fs‖−1∞ , we simply find β∗s = ‖fs‖−1∞ ≥ 1/(2qMψs) and
as a consequence, we obtain
2qMψs
2qµ2∞(ψs)
. (4.16)
A similar relation can be found in the case of a pre-smoothing iteration. Neverthe-
less, since it is enough to prove one between the inf-min conditions, we focus our at-
tention on condition (4.16). So, to enforce the inf–min condition (4.12), it is enough
to prove the existence of an absolute constant L such that µ∞(ψs) 6 L < +∞
uniformly in order to deduce that ‖MGM0‖A0 6
1− (2qL2)−1 < 1.
14 C. Tablino Possio
Proposition 4.6. Under the same assumptions of Lemma 4.5, let us define ψs =
[Φps,q]
s(ψ) for every s ∈ N, where Φp,q is the linear operator defined as in (4.15).
Then, there exists a positive polynomial ψ∞ of degree q such that ψs uniformly
converges to ψ∞, and moreover there exists a positive real number L such that
µ∞(ψs) 6 L for any s ∈ N.
Proof. Due to the periodicity and to the cosine expansions of all the involved
functions, the operator Φq,p in (4.15) can be rewritten as
[Φq,p(ψ)] (x) =
(ϕpψ)
+ (ϕpψ)
. (4.17)
The representation of Φq,p in the Fourier basis (see Proposition 4.8 in [1]) leads
to an operator from Rm(q) to Rm(q), m(q) proper constant depending only on q,
which is identical to the irreducible nonnegative matrix Φ̄q in equation (4.14) of
[1], with q + 1 in place of q.
As a consequence, the claimed thesis follows by referring to the Perron–Frobenius
theorem [15, 30] according to the very same proof technique considered in [1]. �
Lastly, by taking into account all the previous results, we can claim the
optimality of the proposed MGM.
Theorem 4.7. Let f̃0 be a even nonnegative trigonometric polynomial vanishing at
0 with order 2q. Let m0 = m > m1 > . . . > ms > . . . > msmin , ms+1 = ms/2. For
any s = 0, . . . ,msmin−1, let P ss+1 be as in Proposition 4.3 with ps(x) = [1+cos(x)]q,
and let Vs,post = Ims − Ams/‖fs‖∞. If we set Am0 = Cm0(f̃0 + c0χ2πZ) with
c0 = f̃0(w
2 ) and we consider b ∈ Cm0 , then the MGM (standard V-cycle)
converges to the solution of Am0x = b and is optimal (in the sense of Definition
1.1).
Proof. Under the quoted assumptions it holds that f̃0(x) = [1− cos(x)]q ψ0(x)
for some positive polynomial ψ0(x). Therefore, it is enough to observe that the
optimal convergence of MGM as stated in Theorem 2.2 is implied by the inf–min
condition (4.12). Thanks to (4.16), the latter is guaranteed if the quantities µ∞(ψs)
are uniformly bounded and this holds true according to Proposition 4.6. �
Now, we consider the case of a generating function f0 with a unique zero at
x0 = π, this being particularly important in applications since the discretization
of certain integral equations leads to matrices belonging to this class. For instance,
the signal restoration leads to the case of f0(π) = 0, while for the super-resolution
problem and image restoration f0(π, π) = 0 is found [5].
By virtue of Lemma 3.3 we simply have that the generating function f1 related
to the first projected matrix uniquely vanishes at 0, i.e. at the first level the
MGM projects a discretized integral problem, into another which is spectrally and
structurally equivalent to a discretized differential problem.
With respect to the optimal convergence, we have that Theorem 2.2 holds true
with δ = min{δ0, δ̄} since δ results to be a constant and independent of m0.
V-cycle optimal convergence for DCT-III matrices 15
More precisely, δ0 is directly related to the finest level and δ̄ is given by the
inf-min condition of the differential problem obtained at the coarser levels. The
latter constant value has been previously shown, while the former can be proven as
follows: we are dealing with f0(x) = (1+cos(x))
qψ0(x) and according to Definition
4.2 we choose p̃0(x) = p0(x) + d0χ2πZ with p0(x) = (1 + cos(x))
q+1 and d0 =
Therefore, an upper bound for the left-hand side R̃(x) in (4.10) is obtained as
R̃(x) ≤
i.e. we can consider γ0 = Mψ0/m
and so that a value δ0 independent of m0 is
found.
5. Numerical experiments
Hereafter, we give numerical evidence of the convergence properties claimed in
the previous sections, both in the case of proposed TGM and MGM (standard
V-cycle), for two types of DCT-III systems with generating functions having zero
at 0 (differential like problems) and at π (integral like problems).
The projectors P ss+1 are chosen as described in §3 in §4 and the Richardson smooth-
ing iterations are used twice in each iteration with ω = 2/‖f‖ and ω = 1/‖f‖
respectively. The iterative procedure is performed until the Euclidean norm of the
relative residual at dimension m0 is greater than 10
−7. Moreover, in the V-cycle,
the exact solution of the system is found by a direct solver when the coarse grid
dimension equals to 16 (162 in the additional two-level tests).
5.1. Case x0 = 0 (differential like problems)
First, we consider the case Am = Cm(f0) with f0(x) = [2 − 2 cos(x)]q, i.e. with a
unique zero at x0 = 0 of order 2q.
As previously outlined, the matrix Cm(f0) is singular, so that the solution of
the rank-1 corrected system is considered, whose matrix is given by Cm(f0) +
(f0(π/m)/m)ee
T , with e = [1, . . . , 1]T . Since the position of the zero x0 = 0 at the
coarser levels is never shifted, then the function ps(x) = [2− 2 cos(π − x)]r in the
projectors is the same at all the subsequent levels s.
To test TGM/MGM linear convergence with rate independent of the size m0 we
tried for different r: according to (3.6a), we must choose r at least equal to 1 if
q = 1 and at least equal to 2 if q = 2, 3, while according to (4.3a) we must always
choose r equal to q. The results are reported in Table 1.
By using tensor arguments, the previous results plainly extend to the multilevel
case. In Table 2 we consider the case of generating function f0(x, y) = f0(x) +
f0(y), that arises in the uniform finite difference discretization of elliptic constant
coefficient differential equations on a square with Neumann boundary conditions,
see e.g [24].
16 C. Tablino Possio
Table 1. Twogrid/Multigrid - 1D Case: f0(x) = [2 − 2 cos(x)]q
and p(x) = [2− 2 cos(π − x)]r .
q = 1 q = 2 q = 3
m0 r = 1 r = 1 r = 2 r = 2 r = 3
16 7 15 13 34 32
32 7 16 15 35 34
64 7 16 16 35 35
128 7 16 16 35 35
256 7 16 16 35 35
512 7 16 16 35 35
q = 1 q = 2 q = 3
m0 r = 1 r = 1 r = 2 r = 2 r = 3
16 1 1 1 1 1
32 7 16 15 34 32
64 7 17 16 35 34
128 7 18 16 35 35
256 7 18 16 35 35
512 7 18 16 35 35
Table 2. Twogrid/Multigrid - 2D Case: f0(x, y) = [2 −
2 cos(x)]q + [2 − 2 cos(y)]q and p(x, y) = [2 − 2 cos(π − x)]r +
[2− 2 cos(π − y)]r.
q = 1 q = 2 q = 3
m0 r = 1 r = 1 r = 2 r = 2 r = 3
162 15 34 30 - -
322 16 36 35 71 67
642 16 36 36 74 73
1282 16 36 36 74 73
2562 16 36 36 74 73
5122 16 36 36 74 73
q = 1 q = 2 q = 3
m0 r = 1 r = 1 r = 2 r = 2 r = 3
162 1 1 1 1 1
322 16 36 35 71 67
642 16 36 36 74 73
1282 16 36 36 74 73
2562 16 37 36 74 73
5122 16 37 36 74 73
5.2. Case x0 = π (integral like problems)
DCT III matrices Am0 = Cm0(f0) whose generating function shows a unique zero
at x0 = π are encountered in solving integral equations, for instance in image
restoration problems with Neumann (reflecting) boundary conditions [18].
According to Lemma 3.3, if x0 = π, then the generating function f1 of the coarser
matrix Am1 = Cm1(f1), m1 = m0/2 has a unique zero at 0, whose order equals
the order of x0 = π with respect to f0 plus two.
It is worth stressing that in such a case the projector at the first level is singular so
that its rank-1 Strang correction is considered. This choice gives rise in a natural
way to the rank-1 correction considered in §5.1. Moreover, starting from the second
coarser level, the new location of the zero is never shifted from 0.
In Table 3 are reported the numerical results both in the unilevel and two-level
case.
6. Computational costs and conclusions
Some remarks about the computational costs are required in order to highlight
the optimality of the proposed procedure.
V-cycle optimal convergence for DCT-III matrices 17
Table 3. Twogrid/Multigrid - 1D Case: f0(x) = 2+2 cos(x) and
p0(x) = 2 − 2 cos(π − x) and 2D Case: f0(x, y) = 4 + 2 cos(x) +
2 cos(y) and p0(x, y) = 4− 2 cos(π − x)− 2 cos(π − y).
1D TGM MGM
16 15 1
32 14 14
64 12 13
128 11 13
256 10 12
512 8 10
2D TGM MGM
162 7 1
322 7 7
642 7 7
1282 7 6
2562 7 6
5122 7 6
Since the matrix Cms(p) appearing in the definition of P
s+1 is banded, the cost
of a matrix vector product involving P ss+1 is O(ms). Therefore, the first condition
in Definition 1.1 is satisfied. In addition, notice that the matrices at every level
(except for the coarsest) are never formed since we need only to store the O(1)
nonzero Fourier coefficients of the related generating function at each level for
matrix-vector multiplications. Thus, the memory requirements are also very low.
With respect to the second condition in Definition 1.1 we stress that the repre-
sentation of Ams+1 = Cms+1(fs+1) can be obtained formally in O(1) operations
by virtue of (3.5). In addition, the roots of fs+1 and their orders are obtained
according to Lemma 3.3 by knowing the roots of fs and their orders. Finally, each
iteration of TGM costs O(m0) operations as Am0 is banded. In conclusion, each
iteration of the proposed TGM requires O(m0) operations.
With regard to MGM, optimality is reached since we have proven that there exists
δ is independent from both m and smin so that the number of required iterations
results uniformly bounded by a constant irrespective of the problem size. In addi-
tion, since each iteration has a computational cost proportional to matrix-vector
product, Definition 1.1 states that such a kind of MGM is optimal.
As a conclusion, we observe that the reported numerical tests in §5 show that
the requirements on the order of zero in the projector could be weakened. Future
works will deals with this topic and with the extension of the convergence analysis
in the case of a general location of the zeros of the generating function.
References
[1] A. Aricò, M. Donatelli, S. Serra-Capizzano, V-cycle optimal convergence for certain
(multilevel) structured linear systems. SIAM J. Matrix Anal. Appl. 26 (2004), no. 1,
186–214.
[2] A. Aricò, M. Donatelli, A V-cycle Multigrid for multilevel matrix algebras: proof of
optimality. Numer. Math. 105 (2007), no. 4, 511–547 (DOI 10.1007/s00211-006-0049-
18 C. Tablino Possio
[3] O. Axelsson , M. Neytcheva, The algebraic multilevel iteration methods—theory and
applications. In Proceedings of the Second International Colloquium on Numerical
Analysis (Plovdiv, 1993), 13–23, VSP, 1994.
[4] R.H. Chan, T.F. Chan, C. Wong, Cosine transform based preconditioners for total
variation minimization problems in image processing. In Iterative Methods in Linear
Algebra, II, V3, IMACS Series in Computational and Applied Mathematics, Proceed-
ings of the Second IMACS International Symposium on Iterative Methods in Linear
Algebra, Bulgaria, 1995, 311–329.
[5] R.H. Chan, M. Donatelli, S. Serra-Capizzano, C. Tablino-Possio, Application of multi-
grid techniques to image restoration problems. In Proceedings of SPIE - Session: Ad-
vanced Signal Processing: Algorithms, Architectures, and Implementations XII, Vol.
4791 (2002), F. Luk Ed., 210-221.
[6] R.H. Chan and M.K. Ng, Conjugate gradient methods for Toeplitz systems. SIAM
Rev. 38 (1996), no. 3, 427–482.
[7] R.H. Chan, S. Serra-Capizzano, C. Tablino-Possio, Two-grid methods for banded linear
systems from DCT III algebra. Numer. Linear Algebra Appl. 12 (2005), no. 2-3, 241–
[8] G. Fiorentino, S. Serra-Capizzano, Multigrid methods for Toeplitz matrices. Calcolo
28 (1991), no. 3-4, 283–305.
[9] G. Fiorentino, S. Serra-Capizzano, Multigrid methods for symmetric positive definite
block Toeplitz matrices with nonnegative generating functions. SIAM J. Sci. Comput.
17 (1996), no. 5, 1068–1081.
[10] R. Fischer, T. Huckle, Multigrid methods for anisotropic BTTB systems. Linear
Algebra Appl. 417 (2006), no. 2-3, 314–334.
[11] W. Hackbusch,Multigrid methods and applications. Springer Series in Computational
Mathematics, 4. Springer-Verlag, 1985.
[12] T. Huckle, J. Staudacher, Multigrid preconditioning and Toeplitz matrices. Electron.
Trans. Numer. Anal. 13 (2002), 81–105.
[13] T. Huckle, J. Staudacher, Multigrid methods for block Toeplitz matrices with small
size blocks. BIT 46 (2006), no. 1, 61–83.
[14] X.Q. Jin, Developments and applications of block Toeplitz iterative solvers. Combi-
natorics and Computer Science, 2. Kluwer Academic Publishers Group, Dordrecht;
Science Press, Beijing, 2002.
[15] D.G. Luenberger, Introduction to Dynamic Systems: Theory, Models, and Appli-
cations, John Wiley & Sons Inc., 1979.
[16] M.K. Ng, Iterative methods for Toeplitz systems. Numerical Mathematics and Scien-
tific Computation. Oxford University Press, 2004.
[17] M.K. Ng, R.H. Chan, T.F. Chan, A.M. Yip, Cosine transform preconditioners for
high resolution image reconstruction. Conference Celebrating the 60th Birthday of
Robert J. Plemmons (Winston-Salem, NC, 1999). Linear Algebra Appl. 316 (2000),
no. 1-3, 89–104.
[18] M.K. Ng, R.H. Chan, W.C. Tang, A fast algorithm for deblurring models with Neu-
mann boundary conditions. SIAM J. Sci. Comput. 21 (1999), no. 3, 851–866.
V-cycle optimal convergence for DCT-III matrices 19
[19] M.K. Ng, S. Serra-Capizzano, C. Tablino-Possio. Numerical behaviour of multigrid
methods for symmetric sinc-Galerkin systems. Numer. Linear Algebra Appl. 12 (2005),
no. 2-3, 261–269.
[20] D. Noutsos, S. Serra-Capizzano, P. Vassalos, Matrix algebra preconditioners for mul-
tilevel Toeplitz systems do not insure optimal convergence rate. Theoret. Comput. Sci.
315 (2004), no. 2-3, 557–579.
[21] J. Ruge, K. Stüben, Algebraic multigrid. In Frontiers in Applied Mathematics: Multi-
grid Methods.SIAM, 1987, 73–130.
[22] S. Serra Capizzano, Convergence analysis of two-grid methods for elliptic Toeplitz
and PDEs matrix-sequences. Numer. Math. 92 (2002), no. 3, 433–465.
[23] S. Serra-Capizzano, Matrix algebra preconditioners for multilevel Toeplitz matrices
are not superlinear. Special issue on structured and infinite systems of linear equations.
Linear Algebra Appl. 343-344 (2002), 303–319.
[24] S. Serra Capizzano, C. Tablino Possio, Spectral and structural analysis of high preci-
sion finite difference matrices for elliptic operators. Linear Algebra Appl. 293 (1999),
no. 1-3, 85–131.
[25] S. Serra-Capizzano, C. Tablino-Possio, Multigrid methods for multilevel circulant
matrices. SIAM J. Sci. Comput. 26 (2004), no. 1, 55–85.
[26] S. Serra-Capizzano and E. Tyrtyshnikov, How to prove that a preconditioner cannot
be superlinear. Math. Comp. 72 (2003), no. 243, 1305–1316.
[27] H. Sun, X. Jin, Q. Chang, Convergence of the multigrid method for ill-conditioned
block Toeplitz systems. BIT 41 (2001), no. 1, 179–190.
[28] U. Trottenberg, C.W. Oosterlee and A. Schüller, Multigrid. With contributions by
A. Brandt, P. Oswald and K. Stb̈en. Academic Press, Inc., 2001.
[29] E. Tyrtyshnikov, Circulant preconditioners with unbounded inverses. Linear Algebra
Appl. 216 (1995), 1–23.
[30] R.S. Varga, Matrix Iterative Analysis. Prentice-Hall, Inc., Englewood Cliffs, 1962.
C. Tablino Possio
Dipartimento di Matematica e Applicazioni,
Università di Milano Bicocca,
via Cozzi 53
20125 Milano
Italy
e-mail: cristina.tablinopossio@unimib.it
	1. Introduction
	2. Two-grid and Multi-grid methods
	3. Two-grid and Multi-grid methods for DCT III matrices
	4. V-cycle optimal convergence
	5. Numerical experiments
	5.1. Case x0=0 (differential like problems)
	5.2. Case x0= (integral like problems)
	6. Computational costs and conclusions
	References
