E. H. A Y A C H O U R (Lille)
AVOIDING LOOK-AHEAD IN THE LANCZOS METHOD
AND PAD´ E APPROXIMATION
Abstract. In the non-normal case, it is possible to use various look-ahead strategies for computing the elements of a family of regular orthogonal poly- nomials. These strategies consist in jumping over non-existing and singu- lar orthogonal polynomials by solving triangular linear systems. We show how to avoid them by using a new method called ALA (Avoiding Look- Ahead), for which we give three principal implementations. The application of ALA to Pad´e approximation, extrapolation methods and Lanczos method for solving systems of linear equations is discussed.
1. Introduction. A Hankel system comes up implicitly in the Lanczos method, in Pad´e approximation and in extrapolation methods. The princi- pal submatrices of a Hankel matrix are Hankel matrices of linear systems which are solved by using orthogonal polynomials. It is well known that these orthogonal polynomials satisfy a three-term recurrence relation. When some of them are singular, a breakdown (or a so-called true breakdown [8]) problem occurs in this recurrence relation. To avoid such a problem, Draux [16] has shown how to compute regular orthogonal polynomials by using look-ahead strategies. A look-ahead strategy consists in jumping over the non-existing orthogonal polynomials. Draux and Van Ingelandt applied this technique to Pad´e approximation in [17] where they give algorithms which allow moving in the Pad´e table along a diagonal, a row, a staircase consisting of two adjacent diagonals and a sawtooth consisting of two adjacent rows.
Gutknecht and Hochbruck used the Levinson–Schur type recurrences with look-ahead strategies for computing Pad´e approximants [23, 22]. Brezinski, Redivo Zaglia and Sadok have applied these look-ahead strategies to the
1991 Mathematics Subject Classification: 42C05, 41A21, 65F10, 65B05, 65B10.
Key words and phrases : orthogonal and biorthogonal polynomials, Lanczos method, Pad´e approximation, extrapolation methods.
[33]
Lanczos method [10, 12]. These strategies have also been applied to QMR by Freund, Gutknecht and Nachtigal [18, 25].
In this paper, we will show how to substitute new intermediate polyno- mials instead of look-ahead strategies. These intermediate polynomials are biorthogonal and satisfy a simple three-term recurrence relation. In addi- tion, they can be considered as an alternative for orthogonal polynomials which are singular or non-existent.
For a given integer n ∈ Z, we consider the linear functional C
(n)defined on the space C[X] of polynomials by C
(n)(x
i) = c
n+i. By convention, we set c
i= 0 when i < 0. We denote by H
kθnthe following determinant:
H
kθn=
c
θn(0)+nc
θn(0)+n+1. . . c
θn(0)+n+k−1c
θn(1)+nc
θn(1)+n+1. . . c
θn(1)+n+k−1.. . .. . .. .
c
θn(k−1)+nc
θn(k−1)+n+1. . . c
θn(k−1)+n+k−1where θ
nis a permutation of N recursively defined by associating with every j ∈ N the smallest integer θ
n(j) satisfying H
j+1θn6= 0. So, θ
n(0) = i
0if i
0is the smallest integer such that c
i0+n6= 0, θ
n(1) is the smallest integer such that H
2θn6= 0, and so on.
2. Orthogonality. For a fixed integer n, let {P
i(n)}
ibe the family of orthogonal polynomials such that, for all i, P
i(n)has degree i and
(1) C
(n)(x
jP
i(n)(x)) = 0 for j = 0, 1, . . . , i − 1.
For every i and n, if the set of all solutions of (1) is a subspace of dimen- sion 1, then P
i(n)is called regular. The explicit expression of each orthogonal polynomial P
i(n)is given in [8]:
P
i(n)(x) =
c
nc
n+1. . . c
n+ic
n+1c
n+2. . . c
n+i+1.. . .. . .. . c
n+i−1c
n+i. . . c
n+2i−11 x . . . x
i/d
(n)i,
where d
(n)iis the determinant
c
nc
n+1. . . c
n+ic
n+1c
n+2. . . c
n+i+1.. . .. . .. . c
n+i−1c
n+i. . . c
n+2i−1a
n,i0a
n,i1. . . a
n,ii.
Each choice of the coefficients {a
n,ij}
j=0,...,icorresponds to a different nor-
malization of the orthogonal polynomial P
i(n). In the sequel, we will examine three normalizations:
1. In Pad´e approximation, we choose a
n,i0= a
n,i1= . . . = a
n,ii−1= 0 and a
n,ii= 1. Thus P
i(n)is a monic polynomial of degree i.
2. For the Lanczos method, we set a
n,i1= a
n,i2= . . . = a
n,ii= 0 and a
n,i0= 1, which is equivalent to the condition P
i(n)(0) = 1.
3. In extrapolation methods, we choose a
n,i0= a
n,i1= . . . = a
n,ii= 1, which corresponds to P
i(n)(1) = 1.
As we can see from the above explicit expression of P
i(n), the determinant d
(n)ican be zero. This depends on the values assigned to a
n,i0, a
n,i1, . . . , a
n,ii, that is, on the normalization. When P
i(n)is singular, d
(n)iis zero. In this situation, we say that there is a breakdown. The aim of this work is to introduce new regular biorthogonal polynomials with some normalization and to use them in the computation of regular orthogonal polynomials in order to avoid breakdown problems.
Let {P
iθn}
n,ibe a family of monic polynomials such that, for all i, P
iθnhas degree i and
(2) C
(n)(x
θn(j)P
iθn) = 0, j = 0, 1, . . . , i − 1.
The family {P
iθn}
n,icontains all the monic regular orthogonal polynomials with respect to C
(n). The explicit expression of each polynomial P
iθnis
P
iθn(x) =
c
θn(0)+nc
θn(0)+n+1. . . c
θn(0)+n+ic
θn(1)+nc
θn(1)+n+1. . . c
θn(1)+n+i.. . .. . .. .
c
θn(i−1)+nc
θn(i−1)+n+1. . . c
θn(i−1)+n+i1 x . . . x
i/H
iθn.
This shows that P
iθnis a regular orthogonal polynomial if and only if θ
n({0, 1, . . . , i − 1}) = {0, 1, . . . , i − 1}.
In particular, when θ
nis the identity, we recover the explicit expressions of adjacent orthogonal polynomials which are studied, in the normal case, in [5]. When P
iθnis not orthogonal with respect to C
(n), it is, in fact, a biorthogonal polynomial, as defined in [2].
3. Recurrence relations. Assume that there exists a regular orthog- onal polynomial P
kθnsuch that
C
(n)(x
kP
kθn) = 0.
From the explicit expression of P
kθn, we get C
(n)(x
θn(k)P
kθn) = H
k+1θn/H
kθn6= 0. In the following theorem, we see how to compute the biorthogonal polynomials P
k+1θn, P
k+2θn, . . . , P
θθnn(k)of the family {P
iθn}
i, for a fixed inte- ger n.
Theorem 3.1. For i ∈ {k, . . . , θ
n(k) − 1}, we have (i) P
i+1θn= xP
iθn+ α
iP
kθnwith
α
i= −C
(n)(x
θn(k)+1P
iθn)/C
(n)(x
θn(k)P
kθn),
(ii) C
(n)(x
jP
i+1θn) = 0 for j = 0, . . . , θ
n(k) and j = θ
n(k) + k − (i + 1), (iii) θ
n(i + 1) = θ
n(k) + k − (i + 1) and
C
(n)(x
θn(i+1)P
i+1θn) = C
(n)(x
θn(k)P
kθn) 6= 0.
P r o o f. The proof is by induction on i from i = k. It consists in proving that xP
iθn+ α
iP
kθnsatisfies the orthogonality condition (2) for P
i+1θn.
This theorem shows that P
θ(k)+1θnis a regular orthogonal polynomial and that P
kθndivides P
iθnfor i = k, k + 1, . . . , θ
n(k).
Theorem 3.2. Let P
kθn0
be the regular orthogonal polynomial of the high- est degree preceding P
kθn. Then P
θθnn(k)+1
can be computed from the recurrence relation
P
θθnn(k)+1
= xP
θθnn(k)
+
θn(k)
X
i=k+1
α
iP
iθn+ α
kP
kθn+ α
k−1P
kθ0nwith
α
k−1= −C
(n)(x
θn(k)P
kθn)/C
(n)(x
θn(k0)P
kθn0
) 6= 0,
α
k= −(C
(n)(x
θn(k)+1P
θθnn(k)) + α
k−1C
(n)(x
θn(k)P
kθ0n))/C
(n)(x
θn(k)P
kθn), α
i= C
(n)(x
θn(i)P
kθn0
)/C
(n)(x
θn(k0)P
kθn0
), i = k + 1, . . . , θ(k).
P r o o f. Since P
kθn0
is the regular orthogonal polynomial of the highest degree preceding P
kθn, we have θ
n(k
0) = k − 1. For fixed n, the set {P
iθn}
iis a basis of C[X], so we can write
(3) P
θθnn(k)+1= xP
θθnn(k)+
θn(k)
X
i=0
α
iP
iθn.
Multiplying (3) by x
θn(j)and applying C
(n)gives the expressions for α
iand shows that (3) is equivalent to
(4) P
θθnn(k)+1
= xP
θθnn(k)
+
θn(k)
X
i=k+1
α
iP
iθn+ α
kP
kθn+ α
k0P
kθn0
.
By Theorems 3.1 and 3.2, there exists a polynomial W
θn(k)−k+1(x) of de- gree θ
n(k)−k+1 such that P
θ(k)+1θn(x) = W
θn(k)−k+1(x)P
kθn(x)+α
k−1P
kθn0
(x).
It is sufficient to remark that P
kdivides P
iθnfor i = k, . . . , θ
n(k). Different proofs of this property were given by Draux [16], Gragg and Lindquist [19]
and Gutknecht [21]. Notice that our proof is shorter and simpler than that of Gutknecht [21].
The polynomials P
qθncan be displayed in an array called the table P . We suppose that this table P contains a square block of order θ
n(k) − k at its kth column. This can be illustrated by the scheme
P
kθnP
k+1θn−1. . . P
kθn+1P
k+1θn. . . P
kθn+2P
k+1θn+1. . . .. . .. . . ..
where P
kθnis regular. From the preceding results, we obtain the two relations ( P
kθm+1= P
kθm− e
θkmP
θθm+1m+1(k−1)
,
e
θkm= C
(m)(x
kP
kθm)/C
(m+1)(x
k−1P
θθm+1m+1(k−1)
), for m = n, . . . , n + θ
n(k) − k, and
( P
i+1θm−1= xP
iθm− q
i+1θm−1P
θθm−1m−1(i)
, q
θi+1m−1= C
(m)(x
iP
iθm)/C
(m−1)(x
iP
θθm−1m−1(i)
),
for i = k, . . . , θ
n(k), m = n + k − i. These relations yield some properties of blocks of the table P :
Theorem 3.3 For every n ∈ Z, if the table P contains a block at its kth column as described above, then
P
kθn+i= P
kθn, P
k+iθn−i= x
iP
kθn, i = 0, . . . , θ
n(k) − k,
with θ
n−i(k + i) = θ
n(k) and θ
n+i(k) = θ
n(k) − i for i = 0, . . . , θ
n(k) − k.
This was proved by Draux in [16]. Here, we only made the connection with the permutation θ
n, which simplifies the recurrence relations. We also note that a simple proof of this theorem can be deduced from (5) and (6).
The new biorthogonal polynomials defined above are displayed inside the blocks of the table P .
Theorem 3.4. For every n ∈ Z, if the table P has a block at its kth
column as described above, then for i = k, . . . , θ
n(k) and m = 0, . . . , θ
n(k)− i,
we have
P
iθn+m= P
iθn, P
i+mθn−m= x
mP
iθn, θ
n−m(i + m) = θ
n(i), θ
n+m(i) = θ
n(i) − m.
P r o o f. We use the properties of the permutation θ
ngiven in Theorems 3.3 and 3.1. Indeed, from Theorem 3.1, θ
n−m(i + m) = θ
n−m(i + m − 1) + 1, and from Theorem 3.3, θ
n−m(i + m − 1) + 1 = θ
n−m+1(i + m − 1). Thus, by applying Theorem 3.1 and then Theorem 3.3, m times, we deduce that θ
n−m(i + m) = θ
n(i).
Theorem 3.3 shows that θ
n+m(i) = θ
n+m−1(i) − 1. By applying Theo- rem 3.3, m times, we get θ
n+m(i) = θ
n(i) − m.
These relations between the permutations θ
jshow that P
iθnand x
mP
iθnhave the properties of P
iθn+mand P
i+mθn−mrespectively, so P
iθn+m= P
iθnand P
i+mθn−m= x
mP
iθn.
Theorem 3.4 is a generalization of Theorem 3.3.
By using the recurrence relations of Theorems 3.1 and 3.2, we can derive an algorithm for computing the regular orthogonal polynomials with respect to the functional C
(n). Actually, this algorithm allows one to move along a diagonal of the table P . It makes use of the intermediate biorthogonal polynomials for computing the regular orthogonal ones. The procedure is called ALA (Avoiding Look-Ahead strategy).
3.1. Implementation of ALA. Define a symmetric bilinear form g
1on C [X] by
g
1(ψ, ϕ) = C(ψϕ), ∀ψ, ϕ ∈ C[X].
For simplicity, we write C and θ instead of C
(n)and θ
nsince n is fixed.
Definition 3.1. Let D = {p
0, p
1, . . .} and Q = {q
0, q
1, . . .} be two sets of polynomials. If g
1(p
i, q
j) = 0 for i 6= j, then we say that D and Q are g
1-biorthogonal.
We consider two bases {v
0, v
1, . . .} and {w
0, w
1, . . .} of C[X] such that, for every integer i, the polynomials v
iand w
iare of degree i. We assume that (C[X], g
1) is regular. A subspace L ×L
′of C[X]×C[X] is called regular if the right-orthogonal of L,
L
⊥g1= {z ∈ C[X] | ∀y ∈ L, g
1(y, z) = 0},
does not contain any element of L
′, and the left-orthogonal of L
′, L
′g1⊥= {y ∈ C[X] | ∀z ∈ L
′, g
1(y, z) = 0}, is such that L ∩ L
′g1⊥= {0}.
As (C[X], g
1) is regular, we can choose two permutations σ and θ of N
such that, for every integer i, the subspace V
iσ× W
iθgenerated by {v
σ(0), . . .
. . . , v
σ(i)} × {w
θ(0), . . . , w
θ(i)} is regular. This choice enables us to build two g
1-biorthogonal bases D = {p
0, p
1, . . .} and Q = {q
0, q
1, . . .} of monic polynomials such that, for every integer i, p
i∈ V
iσ\V
i−1σand q
i∈ W
iθ\W
i−1θ. For more details, see [2].
To each pair of bases {v
0, v
1, . . .} and {w
0, w
1, . . .}, there corresponds another pair of bases which are g
1-biorthogonal. We are interested in the choice of {v
0, v
1, . . .} and {w
0, w
1, . . .} which yields all the regular orthog- onal polynomials with respect to the functional C. This means that the corresponding g
1-biorthogonal bases {p
0= P
0θ, p
1= P
1θ, . . .} and {q
0= Q
θθ(0), q
1= Q
θθ(1), . . .} satisfy Q
θi= P
iθ= P
ifor every monic regular orthog- onal polynomial P
iof degree i.
We will introduce three interesting choices which give three different ways for implementing the ALA method. These choices will be called C1, C2 and C3.
We now give the recurrence relations connecting the polynomials of {P
0θ, P
1θ, . . .} and {Q
θθ(0), Q
θθ(1), . . .}, in order to apply them to the Lan- czos method. This is equivalent to substituting y for the variable x of the polynomials Q
θi, in order to have two biorthogonal bases with respect to the bilinear form g
2defined on C[X] × C[Y ] by g
2(x
i, y
j) = C(x
i+j) for i, j ∈ N.
• C1 is obtained by taking for σ the identity permutation and by choosing recursively the polynomials of the bases {v
0, v
1, . . .} and {w
0, w
1, . . .} with v
j= xP
j−1θand w
j= x
jfor j = 1, 2, . . . This choice was already studied in Section 2 and before this subsection.
• C2 consists in setting
v
j= xP
j−1θ, j = 1, 2, . . . ,
w
j= x
j−kQ
θk, j = k + 1, . . . , θ(k) + 1,
k = 0, θ(0) + 1, θ(θ(0) + 1) + 1, . . .
The degrees of the polynomials v
jand w
jare equal to their indices. The definitions of P
jθand Q
θjyield the recurrence relations
Q
θj= x
j−kQ
θk, j = k + 1, . . . , θ(k),
(7) (
α
j−1= C(xP
j−1θQ
θθ(k))/C(P
kθQ
θθ(k)),
P
jθ= xP
j−1θ− α
j−1P
kθ, j = k + 1, . . . , θ(k), (8)
α
θ(k)= C(xP
θ(k)θQ
θθ(k))/C(P
kθQ
θθ(k)), β
θ(k)= C(P
θ(k)θQ
θk)/C(P
θ(k−1)θQ
θk−1), P
θ(k)+1θ= xP
θ(k)θ− α
θ(k)P
kθ− β
θ(k)P
θ(k−1)θ, Q
θθ(k)+1= xQ
θθ(k)−
θ(k)
X
i=k+1
α
θ(i)Q
θi− α
θ(k)Q
θk− β
θ(k)Q
θθ(k−1).
(9)
For every i, if the polynomials Q
θiand P
iθare both orthogonal, then Q
θi= P
iθ. Replacing Q
θiby P
iθif Q
θi= P
iθ, these three equations lead to an implementation of the ALA method where only three vectors need to be stored.
The coefficients α
θ(i)of the relation which gives Q
θθ(k)+1can also be computed by using the polynomials of the set {P
k−1′θ, P
k′θ, . . . , P
θ(k)′θ, P
θ(k)+1′θ} which is g
1-biorthogonal to {xQ
θθ(k), Q
θθ(k), Q
θθ(k)−1, . . . , Q
θk, Q
θθ(k−1)} with P
k−1′θ= P
k−1θ. The computation of the polynomials P
i′θis via the following relation which connects them to P
iθ:
(10)
( λ
i= C(P
iθxQ
θθ(k))/C(P
kθQ
θθ(k)),
P
i′θ= P
iθ− λ
iP
k−1θ, i = k, k + 1, . . . , θ(k) + 1.
Thanks to these polynomials, the expression for α
θ(i)is α
θ(i)= −β
θ(k)C(Q
θθ(k−1)P
θ(i)′θ)/C(Q
θiP
θ(i)′θ).
If we only need to compute P
jθ, then we can use the simpler relation
(11)
α
j= C(x
θ(k)−k+1P
jθP
kθ)/C(x
θ(k)−kP
kθP
kθ), β
j= C(P
jθP
kθ)/C(P
θ(k−1)θP
k−1θ),
P
j+1θ= xP
jθ− α
jP
kθ− β
jP
θ(k−1)θ,
j = k, k + 1, . . . , θ(k),
for k = 0, θ(0) + 1, θ(θ(0) + 1) + 1, . . . The initializations of this recurrence relation are P
0θ= 1 and P
−1θ= 0 with θ(−1) = −1.
• C3 consists in taking
v
j= x
j−kP
kθ, j = k + 1, k + 2, . . . , n
k, v
j+1= xP
jθ, j = n
k, n
k+ 1, . . . , θ(k), w
j= x
j−kQ
θk, j = k + 1, k + 2, . . . , n
k, w
j+1= xQ
θj, j = n
k, n
k+ 1, . . . , θ(k),
for k = 0, θ(0) + 1, θ(θ(0) + 1) + 1, . . . , with n
k= ⌊(θ(k) + k + 1)/2⌋. The degrees of v
jand w
jare equal to their indices. For a complete study of this choice, see [2].
For C2, we deduce from (11) the following theorem which generalizes the classical recurrence relation for regular orthogonal polynomials.
Theorem 3.5. Every regular orthogonal polynomial P
θ(k)+1θsatisfies a recurrence relation of the form
P
θ(k)+1θ= xP
θ(k)θ− α
θ(k)P
kθ− β
θ(k)P
θ(k−1)θ,
where P
kθis the regular orthogonal polynomial of the highest degree preceding
P
θ(k)+1θ. The degrees of P
θ(k)+1θ, P
θ(k)θ, P
kθand P
θ(k−1)θare equal to their lower indices.
In the following, we are most interested in the application of C2 because of its particular characteristics which are detailed in [2].
4. Application to the Lanczos method. Let us begin by describing the Lanczos method following [13, 14].
4.1. Description. We want to find the solution of the linear system Ax = b, where A ∈ C
n×nis supposed to be non-singular, b ∈ C
nand x ∈ C
n.
Let x
0and y
0be arbitrary vectors in C
nand define two sequences (x
k)
kand (r
k)
kof vectors by
x
k− x
0∈ K
k(A, r
0), (12)
r
k= b − Ax
k⊥ K
k(A
∗, y
0), (13)
where K
k(A, r) = span(r, Ar, . . . , A
k−1r) and A
∗denotes the conjugate transpose of the matrix A.
The Lanczos method is completely defined by (12) and (13). It con- sists recursively in projecting the initial residual r
0on the Krylov space K
k(A, Ar
0), orthogonally to K
k(A
∗, y
0) with respect to the Hermitian prod- uct h·, ·i of C
n. Here h·, ·i replaces the form g
1introduced before. From (12), we can write
(14) x
k− x
0= −α
1r
0− α
2Ar
0− . . . − α
kA
k−1r
0. Multiplying (14) by A and subtracting b, we obtain
r
k= r
0+ α
1Ar
0+ . . . + α
kA
kr
0. (13) implies
(15) hr
k, A
∗iy
0i = 0 for i = 0, . . . , k − 1.
If we consider the polynomial P
k(ξ) = 1 + α
1ξ + . . . + α
kξ
k, then r
k= P
k(A)r
0. Let us now define the linear functional C on C[X] by C(ξ
i) = c
i= hA
ir
0, y
0i, i = 0, 1, . . . , and the functional C
(1)by C
(1)(ξ
i) = C(ξ
i+1), i = 0, 1, . . . The polynomial P
ksatisfies
C(ξ
iP
k(ξ)) = 0 for i = 0, . . . , k − 1, P
k(0) = 1.
So, P
kis a formal orthogonal polynomial with respect to the linear functional C normalized by the condition P
k(0) = 1.
Let P
k(1)be the monic polynomial of degree k satisfying C
(1)(ξ
iP
k(1)(ξ)) = 0 for i = 0, . . . , k − 1.
(P
k(1))
kand (P
k)
kare called adjacent families [29]. We can easily see that,
for each k ∈N
∗, P
kand P
k(1)exist and are unique if and only if the Hankel
determinant
H
k(1)=
c
1c
2. . . c
kc
2c
3. . . c
k+1.. . .. . .. . c
kc
k+1. . . c
2k−1is different from zero. In order to define uniquely the two sequences (P
kθ)
kand (P
kθ1)
kwith only one permutation θ, (P
kθ)
kand (P
kθ1)
kwill be normal- ized by P
kθ(0) = 1 and P
kθ1monic of degree k.
Even if P
kθis not orthogonal, we set r
k= P
kθ(A)r
0. The polynomial P
kθsatisfies
C(ξ
θ(i)P
kθ(ξ)) = 0 for i = 0, . . . , k − 1, P
kθ(0) = 1.
Consequently, (15) becomes
(16) hr
k, A
∗θ(i)y
0i = 0 for i = 0, . . . , k − 1.
(16) is equivalent to the linear system
(S)
α
1c
θ(0)+1+ α
2c
θ(0)+2+ . . . + α
kc
θ(0)+k= −c
θ(0), α
1c
θ(1)+1+ α
2c
θ(1)+2+ . . . + α
kc
θ(1)+k= −c
θ(1), . . .
α
1c
θ(k−1)+1+ α
2c
θ(k−1)+2+ . . . + α
kc
θ(k−1)+k= −c
θ(k−1).
According to the definition of θ, the determinant of (S) is H
kθ16= 0. So, (S) has a unique solution.
A survey of the various algorithms for implementing the Lanczos method is given in [14]. Here, we only present the application of ALA to Lan- czos/Orthodir which is described in [24, 31].
4.2. Lanczos/Orthodir. Several Lanczos/Orthodir type algorithms were given, for example, in [14]. In particular, we cite the algorithm known as Biodir [21].
According to C2, we apply ALA to Lanczos/Orthodir. This is also equiv- alent to applying ALA to Biodir.
In order to compute P
k+1, we use the formula (17)
( P
i+1θ= P
iθ− λ
iξP
iθ1,
λ
i= C(P
kθQ
θθ(i)1)/C(ξP
kθ1Q
θθ(k)1), i = k, k + 1, . . . , θ(k).
Its proof is by induction from i = k to i = θ(k), it consists in proving that
P
iθ− λ
iξP
iθ1satisfies the orthogonality condition (2) for P
i+1θ. This formula
requires the knowledge of the polynomials P
iθ1, Q
θi1. These polynomials are
obtained from (7)–(9) if we substitute C
(1)for C and θ
1for θ. Therefore,
to apply ALA to Lanczos/Orthodir according to C2, we use the formulas
(7)–(9) and (17).
Now, we are able to give an algorithm which allows avoiding the look- ahead strategy in Biodir. This algorithm consists of three steps:
• initialization,
• determination of the next existing regular orthogonal polynomial, which is equivalent to determining σ(k) at the iteration k,
• computation of the iterates x
k+1, z
k+1and the residual vector r
k+1at the iteration k.
Set z
k= P
kθ1(A)r
0and y
k= Q
θk1(A
∗)y
0for k = 0, 1, . . . Algorithm 1
• Step 1 (Initialization): Choose x
0and y
0arbitrary in C
n, set r
0= b − Ax
0, z
0= r
0, z
−1= y
−1= (0, 0, . . . , 0)
t, h
−1= 1, θ(−1) = −1 and k = 0.
• Step 2 (the determination of σ(k)):
1 i = 0
2 e
i= hy
k+i, r
ki y
k+i+1= A
∗y
k+ih
k+i= hy
k+i+1, z
ki
if |h
k+i| < tol for some tolerance tol, then i = i + 1, go to 2
end (if) θ(k) = k + i.
• Step 3:
b
k= h
θ(k)/h
k−1for i = k, . . . , θ(k) λ
i= e
θ(k)−i/h
θ(k)x
i+1= x
i+ λ
iz
ir
i+1= r
i− λ
iAz
iβ
i= hy
θ(k)+1, Az
ii/h
θ(k)z
i+1= Az
i− β
iz
ky
θ(k)+1← y
θ(k)+1− β
iy
θ(i)end (for)
z
θ(k)+1← z
θ(k)+1− b
kz
θ(k−1)y
θ(k)+1← y
θ(k)+1− b
ky
θ(k−1)k = θ(k) + 1
go to 1
end.
It is important to notice that, for each iteration of this algorithm, we have a product of A and A
∗by a vector and three inner products. The coding of this algorithm needs the storage of 9 + m vectors, where m = max
k(θ(k) − k + 1) < n.
4.3. Numerical results. First, let us mention that the computations were performed on a computer working with 16 decimal digits in double precision and our tests were run using FORTRAN 77.
Let kr
kk be the residual norm obtained, at iteration k, by Algorithm 1.
The algorithm is stopped at the kth iteration if kr
kk < eps, where eps is a given tolerance.
Example 1 . Consider the example of [12]:
0 0 0 . . . 0 −1 1 0 0 . . . 0 0 0 1 0 . . . 0 0 .. . .. . .. . .. . .. . 0 0 0 . . . 1 0
1 2 3 .. . n
=
−n 1 2 .. . n − 1
.
We take n = 1000 and choose y
0= (1, 0, 0, . . . , 0, 0, 1)
t, x
0= (0, 0, . . . , 0)
t. For tol = 10
−1, 10
−2, . . . , 10
−16, eps = 10
−12, we get
θ(0) = 0, θ(k) = 999 − k for k = 1, . . . , 998, θ(999) = 999.
There is stagnation from k = 1 until iteration k = 998. At the end of this stagnation, we obtain kr
999k = 1.58 · 10
4and kr
1000k = 9.55 · 10
−6.
Example 2 . We consider a matrix obtained from discretization of the elliptic partial differential equation
Lu = f on [0, 1] × [0, 1], where
Lu = −∆u + s ∂u
∂x ,
with Dirichlet boundary conditions u = 0, using a five-point centered finite difference scheme on a uniform 20 × 20 grid with mesh size h = 1/21.
This yields a sparse non-symmetric matrix of order n = 400 with 1920
non-zero elements. We choose s = 10
4. By applying Algorithm 1 to this
matrix with tol = 10
−1, 10
−2, . . . , 10
−16, eps = 10
−8, y
0= (0, 0, . . . , 0, 0, 1)
t,
x
0= (0, 0, . . . , 0)
t, b = (1, 0, 0, . . . , 0)
t, we get
0 10 20 30 40 50 60 70
−9
−8
−7
−6
−5
−4
−3
−2
−1 0
Iterations
log10 of residual norm
figure 1, example 2, n=400
As for the first example, there is stagnation at the beginning. Afterwards, we obtain a good convergence to the exact solution. We also remark that the convergence curve presents some peaks. It is well known that these peaks are characteristic of Lanczos type methods.
Example 3. We consider a matrix arising from discretization of the 3-dimensional partial differential equation
Lu = f on [0, 1] × [0, 1] × [0, 1], where
Lu = −∆u + x ∂u
∂x + y ∂u
∂y + z ∂u
∂z − u,
with Dirichlet boundary conditions u = 0. The operator was discretized
using a seven-point centered finite difference scheme on a uniform 5 × 5 × 5
grid with mesh size h equal to 1/6. This yields a sparse non-symmetric
matrix of order n = 125, with 725 non-zero elements. By using Algorithm 1
with tol = 10
−1, 10
−2, . . . , 10
−16, eps = 10
−14, y
0= (0, 0, . . . , 0, 0, 1)
t,
x
0= (0, 0, . . . , 0)
t, b = (1, 0, 0, . . . , 0)
t, we obtain the following convergence
curve:
0 5 10 15 20 25 30
−15
−10
−5 0
Iterations
log10 of residual norm
figure 2, example 3, n=125
For the values of the permutation θ, we get θ(k) =
12 − k for k = 0, 1, . . . , 12, k for k = 13, 14, . . .
At the beginning, we have stagnation from k = 0 until k = 12. After this stagnation, the residual norm converges quickly to zero. At iteration k = 28, we obtain kr
28k = 7.52 · 10
−15.
Let us indicate that we have compared our estimation of the residual norm given by Algorithm 1 with the actual one. This comparison shows that both our estimation and the actual residual norm coincide.
4.4. Application to the non-hermitian Lanczos process. Define the sym- metric bilinear form g
1by g
1(u, w) = w
tu for all u, w ∈ C
n. According to C2, we get the following process.
Process 1 . Choose v
1, v
2∈ C
nand set λ
1p
1= v
1, µ
1q
1= v
2, p
0= 0, k = 1 (λ
1and µ
1are chosen such that kp
1k = kq
1k = 1).
Compute 1 i = 0
2 d
k= g
1(p
k, q
k+i)
if |d
k| < tol (for some tolerance tol), then i = i + 1
µ
k+iq
k+i= A
∗q
k+i−1(µ
k+iis chosen such that kq
k+ik = 1) go to 2
end (if)
θ(k + j) = k + i − j, j = 0, 1, . . . , i q
θ(k)+1= A
∗q
θ(k)for j = k, k + 1, . . . , θ(k):
α
j= g
1(Ap
j, q
θ(k))/g
1(p
k, q
θ(k)) p
j+1= Ap
j− α
jp
kβ
j= g
1(q
θ(k)+1, p
j)/g
1(q
θ(j), p
j) q
θ(k)+1← q
θ(k)+1− β
jq
θ(j)if j = θ(k), then
α
′j= g
1(Ap
j, q
k−1)/g
1(p
θ(k−1), q
k−1) p
j+1← p
j+1− α
′jp
θ(k−1)β
j′= g
1(q
θ(k)+1, p
k−1)/g
1(q
θ(k−1), p
k−1) µ
θ(k)+1q
θ(k)+1← q
θ(k)+1− β
j′q
θ(k−1)(µ
θ(k)+1is chosen such that kq
θ(k)+1k = 1) end (if)
λ
j+1p
j+1← p
j+1(λ
j+1is chosen such that kp
j+1k = 1) end (for)
k = θ(k) + 1, go to 1 end.
For solving a linear system, we use a process which allows us to triangu- larize, tridiagonalize or transform the matrix of the system to another one for which we have to find its inverse, as for example the Hessenberg matrix.
Here, for each iteration k of Process 1, we get the following factorization:
A(p
kp
k+1. . . p
θ(k)) = (p
θ(k−1)p
kp
k+1. . . p
θ(k)p
θ(k)+1)
d
′ktH e
k′where d
′kt= α
′θ(k)(0, 0, . . . , 0, 1) ∈ C
θ(k)−k+1, and the matrix e H
′kwith θ(k) − k + 2 rows and θ(k) − k + 1 columns is
α
kα
k+1. . . α
θ(k)λ
k+1λ
k+2. ..
λ
θ(k)+1
.
Let us discuss the stopping criterion for this process. We consider two
Krylov subspaces W
2= K
n(A, v
1) = span(v
1, Av
1, A
2v
1. . .) and W
3=
K
n(A
∗, v
2) = span(v
2, A
∗v
2, A
∗2v
2, . . .). There are two cases to consider:
• The first one corresponds to not having a breakdown at iteration k if k ≤ min{l, l
′} with l = dim W
2and l
′= dim W
3. This means that the subspace (W
2× W
3, g
1) is regular.
• The second case is the situation where there is a serious incurable breakdown. It corresponds to a breakdown occurring at iteration k with k ≤ min{l, l
′} and it means that (W
2×W
3, g
1) is not regular. Consequently, Process 1 cannot be used and the solution is to make another choice of the vectors v
1and v
2in C
n.
Remark 4.1. This process needs the storage of m + 5 vectors of C
n, where m = max
k{θ(k) − k + 1}. m + 5 is smaller than the number 2m + 4 of vectors used in look-ahead strategies. In the regular case, the classical Lanczos process only needs 6 vectors. This number coincides with m + 5, since in the regular case, θ(k) = k, which implies that m = 1 for all k.
We note that the factorization of Process 1 has also been used by Ziegler in [32, 33], where he talks about a special look-ahead strategy.
Remark 4.2. We have shown how to apply C2 to the Lanczos method.
We can do the same for the CGM-type (Conjugate Gradient Multiplied) methods which have been simultaneously introduced by Brezinski [7] and Gutknecht [20], and which are also known under the name of “product-type methods”. The CGM class contains CGS (Conjugate Gradient Squared) due to Sonneveld [27] and Bi-CGSTAB due to Van Der Vorst [28].
5. Application to Pad´ e approximation. Orthogonal polynomials and their associates implicitly come up in the computation of Pad´e approx- imants. Blocks of a non-normal Pad´e table are due to the non-existence and singularity of some orthogonal polynomials. In this section, we give relations between orthogonal, reciprocal, associated and intermediate poly- nomials introduced in [1], and we show how to apply them to the recursive computation of Pad´e approximants.
Let f be a formal power series f (t) = c
0+ c
1t
1+ c
2t
2+ . . . with c
i∈ C for i ∈ N. We look for a rational fraction
R(t) = Q(t)
P (t) = a
0+ a
1t + . . . + a
pt
pb
0+ b
1t + . . . + b
qt
qwhose power series expansion in ascending powers of t agrees with f as far
as possible, which means that f (t) − R(t) = O(t
p+q+1) (t → 0). Such a
rational fraction is called a Pad´e approximant of f and it is denoted by
[p/q]
f(t). Usually these approximants are displayed in a two-dimensional
array called the Pad´e table. Identical Pad´e approximants can only occur in
square blocks in the Pad´e table. If there is no block, we say that the Pad´e
table is normal. Otherwise, it is called non-normal.
For every n ∈ Z, define the linear functional C
(n)on the space of complex polynomials by C
(n)(x
i) = c
n+iwith the convention that c
i= 0 for i < 0.
C
(n)is associated with the formal power series
f
n(t) = c
n+ c
n+1t + c
n+2t
2+ . . .
5.1. The associated polynomials. For every P
qθn, we consider the associ- ated polynomial
Q
θqn(t) = C
(n)P
qθn(x) − P
qθn(t) x − t
where C
(n)acts on x.
Lemma 5.1. If Q
θknis associated with the polynomial P
kθnof degree k, then
Q
θkn(t) = X
m i=0t
iC
(n−i−1)(P
kθn(x))
where C
(n−i−1)acts on x and m = n + k − 1 − θ
n(0). Q
θkn(t) has degree m if m ≥ 0, otherwise Q
θkn(t) = 0.
P r o o f. Q
θkn(t) is equal to C
(n)[(P
kθn(x) − P
kθn(t))/(x − t)]. By using the equality
1/(x − t) = x
−1X
∞ i=0(x
−1t)
i, we prove that
Q
θkn(t) = C
(n)[P
kθn(x) − P
kθn(t)]x
−1n+k−1
X
i=0
(x
−1t)
i. Finally, since c
i= 0 for i < θ(0), we obtain the result of the lemma.
5.2. The reciprocal orthogonal polynomials. We consider the reciprocal series g of t
−θ(0)f defined by t
−θ(0)f (t)g(t) = 1. We set g(t) = P
∞i=0
d
it
iand we define a functional D
(n)on C[X] by D
(n)(x
i) = d
n+ifor i ∈ N. D
(n)is called the reciprocal functional of C
(n). By convention, we set d
i= c
i= 0 if i < 0. Let η
nbe the permutation associated with the functional D
(n); it is called the reciprocal permutation of θ
n. We remark that the definition of D
(0)= D implies η(0) = η
0(0) = 0. We will find later a relation which gives us the permutation η
nfrom θ
n. The complex numbers d
iare obtained from the equations
c
θ(0)d
0= 1, c
θ(0)d
j+ c
θ(0)+1d
j−1+ . . . + c
θ(0)+jd
0= 0 for j = 1, 2, . . . An orthogonal polynomial with respect to D
(n)is called reciprocal. We denote by {R
ηin}
ithe family of all these reciprocal orthogonal polynomials.
They are useful for the recursive computation of numerators of Pad´e ap-
proximants. In the following theorem, we study the connection between the polynomials of the two families {P
iθn}
i,nand {R
ηin}
i,n.
Theorem 5.1. If one of the polynomials P
kθθ(0)+n+1and R
ηn+k−n+1is regular and orthogonal, then so is the other. The same holds for P
n+kθθ(0)−n+1and R
kηn+1. If P
kθθ(0)+n+1and P
n+kθθ(0)−n+1are regular and orthogonal, then
S
n+kη−n+1= d
0P
kθθ(0)+n+1, Q
θn+kθ(0)−n+1= c
θ(0)R
ηkn+1, n = 1, 2, . . . ,
c
θ(0)R
ηn+k−n+1= P
kθθ(0)+n+1X
n i=0c
θ(0)+ix
n−i+ Q
θkθ(0)+n+1,
d
0P
n+kθθ(0)−n+1= R
ηkn+1X
ni=0
d
ix
n−i+ S
kηn+1,
n = 0, 1, . . .
P r o o f. It is sufficient to remark that f
θ(0)is the reciprocal series of g (this means that C
(θ(0))is the reciprocal functional of D) and then use the results of [5, 16]. When θ(0) = 0, the proof given in [5, 4] of the equalities of this theorem is long. It consists in transforming the determinants of the explicit expressions of the orthogonal polynomials. A simple proof is obtained by using only Lemma 5.1 (see the proof of Theorem 5.2).
From this theorem, it is clear that, for a fixed integer n, R
n+kη−n+1, S
n+kη−n+1, P
kθθ(0)+n+1and Q
θkθ(0)+n+1satisfy the same recurrence relations with differ- ent initializations. The same holds for P
n+kθθ(0)−n+1, Q
θn+kθ(0)−n+1, R
ηkn+1and S
kηn+1. If we set, for every n, k ∈ N, N
kηn+2= c
θ(0)R
ηk−nand N
kη−n+1= c
θ(0)R
ηkn+1, then the Pad´e approximant [p/q]
fcan be written as [p/q]
f= N e
pηp−q+1/ e P
θqθ(0)+p−q+1whenever P
qθθ(0)+p−q+1is regular (see [5]).
We deduce from the preceding results that whether or not there are blocks in the Pad´e table, the numerator of each Pad´e approximant can be computed recursively by using the recurrence relations satisfied by the denominators.
Corollary 5.1. The permutations η
nare connected to θ
nby the follow- ing relations:
η
−n+1(i) = θ
θ(0)+n+1(i − n) + n,
θ
θ(0)−n+1(i) = η
n+1(i − n) + n for i ≥ n, n ≥ 1,
θ
θ(0)−n+1(i) = η
−n+1(i) = n − 1 − i for i = 0, 1, . . . , n − 1.
P r o o f. The knowledge of the degrees of all the regular orthogonal
polynomials implies that of η
nand θ
n, see Theorems 3.1 and 3.2. So, from
the definition of the permutations θ
n, η
nand by using Theorem 5.1, we get
the assertion.
The quantities Q
θn+kθ(0)−n+1and P
kθθ(0)+n+1P
ni=0
c
θ(0)+ix
n−i+ Q
θkθ(0)+n+1intervene in the recursive computation of the numerators of the Pad´e ap- proximants. These quantities do not satisfy the equalities of Theorem 5.1 when P
kθθ(0)+n+1and P
n+kθθ(0)−n+1are not orthogonal. For this reason, we give some properties of them below. For every n ∈ Z and k ∈ N, we consider the monic polynomials R
′kηndefined by
(18) D
(n)(R
′kηnt
ηn(j)) + α
n,kd
ηn(j)−ηn(k)= 0, j = 0, . . . , k − 1, where α
n,kis a constant such that the solution R
k′ηnof (18) is monic. The role of these polynomials is to replace the polynomials R
ηknin the equalities of Theorem 5.1. By substituting R
′kηnfor R
ηkn, the results of Theorem 5.1 are true even if P
kθθ(0)+n+1and P
n+kθθ(0)−n+1are not orthogonal. Clearly, the family {R
′kηn}
n,kis built in such a way that it contains all the regular orthogonal polynomials with respect to the functional D
(n). From the definition of R
k′ηn, we can easily see that the condition for their existence and unicity is the same as for the polynomials R
ηkn. Therefore, for every k and n, the polynomial R
′kηnexists, is unique and has degree k.
Theorem 5.2. We have
Q
θn+kθ(0)−n+1= c
θ(0)R
′kηn+1, S
n+k′η−n+1= d
0P
kθθ(0)+n+1, n = 1, 2, . . . ,
c
θ(0)R
′n+kη−n+1= P
kθθ(0)+n+1X
n i=0c
θ(0)+ix
n−i+ Q
θkθ(0)+n+1,
d
0P
n+kθθ(0)−n+1= R
′kηn+1X
ni=0
d
ix
n−i+ S
k′ηn+1,
n = 0, 1, . . .
P r o o f. We want to prove Q
θn+kθ(0)−n+1= c
θ(0)R
′kηn+1. First assume that θ(0) = 0. In this case, thanks to Lemma 5.1, we have, for j = 0, 1, . . . , k − 1,
D
(n+1)[Q
θn+k−n+1(t)t
ηn+1(j)] = D
(n+1)h X
ki=0
t
i+ηn+1(j)C
(−n−i)(P
n+kθ−n+1(x)) i
=
n+k
X
l=0
a
lX
ki=0
d
i+ηn+1(j)+n+1c
−n−i+lwhere D
(n+1)acts on t, C
(−n−i)acts on x and P
n+kθ−n+1(x) = P
n+k l=0a
lx
l. Since θ(0) is zero, Corollary 5.1 implies θ
−n+1(j + n) = η
n+1(j) + n, and we conclude that
D
(n+1)[Q
θn+k−n+1(t)t
ηn+1(j)] =
n+k
X
l=0
a
lX
l−n i=0d
i+θ−n+1(j+n)+1
