On the Quadratic Convergence of the Falk-Langemeyer Method

Ivan Slapničar^{∗} Vjeran Hari^{†}

Abstract

The Falk–Langemeyer method for solving a real definite generalized eigenvalue problem, Ax = λBx, x ≠ 0, is proved to be quadratically convergent under arbitrary cyclic pivot strategy if the eigenvalues of the problem are simple. The term “quadratically convergent” roughly means that the sum of squares of the off–diagonal elements of matrices from the sequence of matrix pairs generated by the method tends to zero quadratically per cycle.

Key words: generalized eigenvalue problem, Jacobi method, quadratic convergence, asymptotic convergence

AMS(MOS) subject classification. 65F15, 65F30

∗Department of Mathematics, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, University of Split, R. Boškovića b.b., YU-58000 Split, Yugoslavia

†Department of Mathematics, University of Zagreb, pp.187, YU-41000 Zagreb, Yugoslavia

1 Introduction

In this paper we study the asymptotic convergence of the method established in 1960 by S. Falk and P. Langemeyer in [2]. Their method solves generalized eigenvalue problem

Ax = λBx ,  x ≠ 0,  (1)

where A and B are real symmetric matrices of order n such that the pair (A, B) is definite. By definition the pair (A, B) is definite if the matrices A and B are hermitian or real symmetric and there exist real constants a and b such that the matrix a A + b B is positive definite.

The Falk–Langemeyer method is the most commonly used Jacobi–type method for solving problem (1). Its advantages over other methods of solving problem (1) are that it applies to problem (1) for the widest class of starting pairs. Although it is not, in general, the fastest method for solving the given problem, in some cases it is the most appropriate. The QR method [11] is usually several times faster, at least on conventional computers, but it solves problem (1) only if matrix B is positive definite (or positive definitizing shift for the pair is known in advance) and if matrix B is well conditioned for Cholesky decomposition. The Falk–Langemeyer method is superior to the QR method in terms of numerical stability if matrix B is badly conditioned for Cholesky decomposition. It is also superior to the QR method if approximate eigenvectors are known, i.e. if the matrices A and B are almost diagonal. This happens in the course of modeling the parameters of a system where a sequence of matrix pairs differing only slightly from each other has to be reduced. This also happens in various subspace iteration techniques (see [11]). Another reason why Jacobi–type methods have attracted attention recently is that they are adaptable for parallel processing (see [12], [10]).

The Jacobi–type method for solving problem (1) recently proposed by Veselić in [15] is somewhere in between previously mentioned methods in both, speed and requirements. Although Veselić’s method works for definite matrix pairs, a linear combination ρA − λB which is reasonably well conditioned for J–symmetric Cholesky decomposition must be known in advance. This method is one of the implicit methods, i.e. it works only on the eigenvectors matrix, and is therefore approximately two times faster than the Falk–Langemeyer method.

The Jacobi–type method considered by Zimmerman in [19] is closely related to the Falk–Langemeyer method (this is briefly described in Section 3) but requires positive definite matrix B. In [19] the convergence of this method is proved under the assumption that the starting matrices are almost diagonal. The same conclusion holds for the Falk–Langemeyer method as we shall show in this paper.

In [4] Hari studied the asymptotic convergence of complex extension of Zimmerman’s method (also for positive definite B). He showed that his method converges quadratically under the cyclic pivot strategies if the eigenvalues of the problem are simple, while in the case of multiple eigenvalues the method can be modified so that the quadratic convergence persists. We are interested only in cyclic pivot strategies since some of them are amenable for parallel processing.

These results, the informal analysis of the convergence properties of the Falk–Langemeyer method performed by Hari in [7], and the numerical investigation suggested that the Falk–Langemeyer method behaves in the similar fashion. In this paper we prove that the Falk–Langemeyer method is quadratically convergent if the eigenvalues of the problem are simple and the pivot strategy is cyclic. The technique of the proof, originally established by the late J. H. Wilkinson in [16] (cf. [6]), is similar to that used in [4] .

Two main problems that had to be solved are that neither of the matrices A and B has to be positive definite and that the transformation matrices are not orthogonal and therefore difficult to estimate. Both problems were solved using the results about almost diagonal definite matrix pairs from [7].

The paper is organized as follows. In Section 2 we state the known results about almost diagonal definite matrix pairs from [7] to the extent necessary for understanding the rest of the paper. In Section 3 we describe the Falk–Langemeyer method, show that it always works for definite matrix pairs (without use of definitizing shifts), and give its algorithm. We also briefly describe Zimmerman method from [19] and [4] and relate it to the Falk–Langemeyer method. Section 4 is the central section of the paper. We first state the known result about the quadratic convergence of Zimmerman method from [4] and show to what extent can this result be applied to the Falk–Langemeyer method. We introduce measure ε̃_{k} which we use for defining and proving quadratic convergence. Then we prove the quadratic convergence of the Falk–Langemeyer method under the assumptions that the eigenvalues of the problem are simple, the pivot strategy is arbitrary cyclic and the matrices A and B are almost diagonal. At the end we show that the quadratic convergence implies the convergence of Falk–Langemeyer method. In Section 5 we give the quadratic convergence results for parallel and serial strategies, briefly explain the possible modification of the Falk–Langemeyer method in case of multiple eigenvalues, and briefly discuss numerical experiments.

Most of the results presented in the paper are part of an M. S. thesis [13] done under the supervision of professor V. Hari.

We would like to thank professor K. Veselić from Fernuniversität Hagen for his helpful suggestions. We would also like to thank both reviewers for their comments which helped us clarify some important parts of the paper.

2 Almost Diagonal Definite Matrix Pairs

Here we consider the structure of almost diagonal definite matrix pair. We first state some properties of definite matrix pairs. Then we introduce chordal metric for measuring distance between eigenvalues of definite matrix pairs. We define the measures for the almost diagonality of the square matrix and of the pair of square matrices. At the end we state an important theorem from [7]. The theorem and its corollary reveal the structure of almost diagonal definite matrix pairs. All results are given for the general case of hermitian matrices even though in the rest of the paper we shall consider only the case of real symmetric matrices.

Definite matrix pair (A, B) has some important properties:

a) There exists a nonsingular matrix F such that

F^{∗}AF = diag (a_{1}, …, a_{n}) = D_{A}
F^{∗}BF = diag (b_{1}, …, b_{n}) = D_{B}.  (2)

The ratios a_{i}/b_{i} , i = 1, …, n, of real numbers a_{i}, b_{i} are the eigenvalues of the pair (A, B) and are unique to the ordering. If [f_{1}, ..., f_{n}] denotes the partition by columns of F, vectors f_{i} , i = 1, ..., n, are the corresponding eigenvectors. Matrices D_{A} and D_{B} are not uniquely determined by the pair (A, B). In the real symmetric case F^{∗} can be changed to F^{T} in the relation (2).

b) The Crawford constant c(A, B),

c(A, B) = inf{| x^{∗}(A + iB)x | ; x ∈ C^{n}, ∥ x ∥= 1}  (3)

is positive. Therefore, A and B share no common nul–subspace and |a_{i}| + |b_{i}| > 0, i = 1, …, n, independently of the choice of F. Note that the choice x = e_{i} (the i-th coordinate vector) in the relation (3) for i = 1, …, n implies

d_{i} = ^{4} \sqrt{(a_{ii})^{2} + (b_{ii})^{2}} > 0 ,  i = 1, …, n ,  (4)

where A = (a_{ij}) and B = (b_{ij}). Hence the matrix

D = diag (\frac{1}{d_{1}} , …, \frac{1}{d_{n}}) ,  (5)

is positive definite. In the real symmetric case for n ≠ 2 only real vectors x can be taken in the relation (3).

c) There exists a real number φ, such that the matrix B_{φ} from the pair (A_{φ}, B_{φ}),

A_{φ} = A cos φ − B sin φ
B_{φ} = A sin φ + B cos φ ,  (6)

is positive definite. The matrices A and B can be simultaneously diagonalized if and only if the same holds for the matrices A_{φ} and B_{φ}.

The proofs of the above properties are simple (see [14]). If some f_{i} is a vector from the nul–subspace of B, the eigenvalue λ_{i} is infinite. Such eigenvalues are not badly posed because they are zero eigenvalues of the pair (B, A) counting their multiplicities. Hence, it is better to define eigenvalues as pairs of numbers λ_{i} = [a_{i}, b_{i}], i = 1, …, n. It is also necessary to choose a finite metric for measuring the distance between eigenvalues. Such is the chordal metric.

Let R^{2} = R × R and R_{0}^{2} = R^{2} \ {[0, 0]}, where R is the set of real numbers. We say that the pairs [a, b], [c, d] ∈ R_{0}^{2} are equivalent if ad − bc = 0 and write [a, b] ρ [c, d]. It is easily seen that ρ is an equivalence relation on R_{0}^{2}. Let R_{0}^{2} |_{ρ} be the set of equivalence classes. Let λ, μ ∈ R_{0}^{2} |_{ρ} and let [a, b], [c, d] be their representatives, respectively. Chordal distance between [a, b] and [c, d] is defined with the formula

χ([a, b], [c, d]) = \frac{| ad − bc |}{\sqrt{a^{2} + b^{2}} \sqrt{c^{2} + d^{2}}} .

It is easily seen that χ is constant when [a, b] and [c, d] vary over λ and μ, respectively. This defines metric χ̃ : R_{0}^{2} |_{ρ} ×R_{0}^{2} |_{ρ}→ R by χ̃(λ, μ) = χ([a, b], [c, d]) where [a, b] and [c, d] are any representatives of λ and μ, respectively. However, for the sake of simplicity we shall use χ for the both functions χ and χ̃. We see that χ(λ, μ) ≤ 1 for all λ, μ ∈ R_{0}^{2} |_{ρ}. The proof of these and some other properties of the chordal metric can be found in [11],

[14] and [13].

From now on, let n denote the order of the matrices A and B and let p denote the number of distinct eigenvalues of the pair (A, B). We assume that

n ≥ 3 , p ≥ 2 ,

and that the pair (A, B) is definite. Note that if p = 1 then A = λB, so λ_{1} = λ_{2} = ⋯ = λ_{n} = λ and all vectors are eigenvectors of the pair (A, B).

The off–norm of the square matrix A is the quantity

S(A) = \sqrt{∑_{i≠j^{i,j=1}}^{n} | a_{ij} |^{2}} = ∥ A − diag (A) ∥ ,

where ∥ · ∥ denotes the Euclidean matrix norm.

The off–norm of the pair (A, B) is the quantity

ε(A, B) = \sqrt{S^{2}(A) + S^{2}(B)} .  (7)

Where no misunderstanding can arise, ε shall be used instead of ε(A, B). Let

λ_{1} = … = λ_{t_{1}} , λ_{t_{1}+1} = … = λ_{t_{2}} , …, λ_{t_{p−1}+1} = … = λ_{t_{p}} ,  (8)

where

λ_{t_{i}} = [s_{i}, c_{i}] , s_{i}^{2} + c_{i}^{2} = 1 , i = 1, …, p , (9)

be all eigenvalues of the pair (A, B). Thus, we assume that the pair (A, B) has p distinct eigenvalues λ_{t_{1}}, …, λ_{t_{p}} with the appropriate multiplicities

n_{i} = t_{i} − t_{i−1} ,  i = 1, …, p ,  t_{0} = 0 ,  (10)

and the representatives which behave as sine and cosine are chosen. Since p > 2 we can define quantities

δ_{i} = \frac{1}{3} min_{j≠i^{1≤j≤p}} χ(λ_{t_{i}}, λ_{t_{j}}) ,  δ = min_{1≤i≤p} δ_{i} .  (11)

Note that δ > 0.

In the analysis we shall need matrices Ã and B̃ defined as

Ã = DAD ,  B̃ = DBD ,  (12)

where matrix D is defined with relations (5) and (4). Since D is positive definite, the matrices Ã and B̃ are congruent to the matrices A and B, respectively. Let us partition the matrices Ã and B̃,

Ã = [ \begin{matrix}Ã_{11} & ⋯ & Ã_{1p} \\ ⋮ & & ⋮ \\ Ã_{p1} & ⋯ & Ã_{pp} \end{matrix} ] ,  B̃ = [ \begin{matrix}B̃_{11} & ⋯ & B̃_{1p} \\ ⋮ & & ⋮ \\ B̃_{p1} & ⋯ & B̃_{pp} \end{matrix} ] ,  (13)

where Ã_{ii} and B̃_{ii} are diagonal blocks of order n_{i} , i = 1, …, p, and n_{i}’s are defined with relation (10). The relation (13) shall be written as Ã = (Ã_{ij}) and B̃ = (B̃_{ij}).

Let the matrices A and B be partitioned according to the relation (13). Departure from the block–diagonal form of the pair (A, B) is the quantity

τ(A, B) = \sqrt{τ^{2}(A) + τ^{2}(B)} ,

where

τ^{2}(A) = ∑_{i=1}^{p} ∑_{j≠i^{j=1}}^{p} ∥ A_{ij} ∥^{2} , τ^{2}(B) = ∑_{i=1}^{p} ∑_{j≠i^{j=1}}^{p} ∥ B_{ij} ∥^{2} .

Theorem 1 Let (A, B) be a definite pair and let the matrices Ã and B̃ be defined by the relations (12), (5) and (4). If

ε(Ã, B̃) < δ ,  (14)

then there exists a permutation matrix P such that for matrices Ã′ = P^{T} ÃP and B̃′ = P^{T} B̃P, partitioned according to the relation (13), holds

∥ c_{i}Ã′_{ii} − s_{i}B̃′_{ii} ∥ ≤ \frac{1}{δ_{i}} ∑_{j≠i^{j=1}}^{p} ∥ c_{i}Ã′_{ij} − s_{i}B̃′_{ij} ∥^{2} , i = 1, …, p .  (15)

On the both sides of the inequalities (15) the Euclidean matrix norm can be substituted with the spectral norm.

Proof: The proof of this theorem is found in [7]. Q.E.D.

Corollary 2 Let the relation (14) hold for the definite pair (A, B). Then there exists a permutation matrix P such that for the matrices Ã′ = P^{T} ÃP = (ã′_{ij}) , B̃′ = P^{T} B̃P = (b̃′_{ij}), A′ = P^{T} AP = (a′_{ij}) and B′ = P^{T} BP = (b′_{ij}), partitioned according to the relation (13), holds

∑_{i=1}^{p} ∥ c_{i}Ã′_{ii} − s_{i}B̃′_{ii} ∥^{2} ≤ \frac{τ^{4}(Ã′, B̃′)}{2δ^{2}} ,  (16)
∑_{i=1}^{p} ∑_{j=t_{i−1}+1}^{t_{i}} χ^{2}([s_{i}, c_{i}], [a′_{jj}, b′_{jj}]) = ∑_{i=1}^{p} ∑_{j=t_{i−1}+1}^{t_{i}} | c_{i}ã′_{jj} − s_{i}b̃′_{jj} |^{2}

≤ \frac{τ^{4}(Ã′, B̃′)}{2δ^{2}} ,  (17)
χ([s_{i}, c_{i}], [a′_{jj}, b′_{jj}]) = | c_{i}ã′_{jj} − s_{i}b̃′_{jj} | ≤ \frac{τ^{2}(Ã′, B̃′)}{2δ} ,
j = t_{i−1} + 1, …, t_{i} , i = 1, …, p . (18)

Proof: By Theorem 1 there exists a permutation matrix P such that the relation (15) holds for the matrices Ã′ and B̃′. The Cauchy–Schwarz inequality implies

∥ c_{i}Ã′_{ij} − s_{i}B̃′_{ij} ∥^{2} ≤ (| c_{i} |∥ Ã′_{ij} ∥ + | s_{i} |∥ B̃′_{ij} ∥)^{2}
≤ ∥ Ã′_{ij} ∥^{2} + ∥ B̃′_{ij} ∥^{2} ,  i ≠ j .  (19)

From the relations (15) and (19), the definition of τ(Ã′, B̃′) , and the symmetry of matrices Ã′ and B̃′ follows

∥ c_{i}Ã′_{ii} − s_{i}B̃′_{ii} ∥ ≤ \frac{1}{δ_{i}} ∑_{j≠i^{j=1}}^{p}(∥ Ã′_{ij} ∥^{2} + ∥ B̃′_{ij} ∥^{2})

≤ \frac{1}{2δ} τ^{2}(Ã′, B̃′) ,  i = 1, …, p .  (20)

Finally, the relations (15), (19) and (20) and the definition of τ(Ã′, B̃′) imply

∑_{i=1}^{p} ∥ c_{i}Ã′_{ii} − s_{i}B̃′_{ii} ∥^{2} ≤ \frac{1}{2δ} τ^{2}(Ã′, B̃′) ∑_{i=1}^{p} ∑_{j≠i^{j=1}}^{p} \frac{1}{δ_{i}} (∥ Ã′_{ij} ∥^{2} + ∥ B̃′_{ij} ∥^{2})

≤ \frac{1}{2δ^{2}} τ^{4}(Ã′, B̃′) ,

which completes the proof of the relation (16).

The equalities in the relations (17) and (18) follow from the definition of the chordal metric and the fact that it does not depend upon the choice of the representatives. Inequality in the relation (17) now follows from the relation (16) and inequality in the relation (18) from the relation (20). Q.E.D.

Theorem 1 and Corollary 2 reveal the structure of almost diagonal definite matrix pairs in both the hermitian and the real symmetric case. The relation (18) implies that for i = 1, …, p, pairs [a′_{jj}, b′_{jj}] , j ∈ {t_{i−1} + 1, …, t_{i}} approximate the eigenvalues λ_{t_{i}} with an error of order of magnitude τ^{2}(Ã′, B̃′) in the chordal metric. The relation (16) implies that the blocks Ã′_{ii} and B̃′_{ii}, i = 1, …, p, are proportional with the proportionality constants being λ_{t_{i}} also with the error of order τ^{2}(Ã′, B̃′). This proportionality becomes apparent when τ(Ã′, B̃′) is small enough compared to δ. Note that the relations (15) and (18) do not imply that the off-diagonal elements of blocks Ã′_{ii} and B̃′_{ii} tend to zero together with τ(Ã′, B̃′). The relation (15) shows that for fixed i the proportionality of the blocks Ã′_{ii} and B̃′_{ii} depends on the local separation δ_{i} of the eigenvalue λ_{t_{i}} from other eigenvalues and on quantities ∥ c_{i}Ã′_{ij} − s_{i}B̃′_{ij} ∥^{2} , j = 1, …, p , j ≠ i.

æ

3 The Falk–Langemeyer method

In this section we define the Falk–Langemeyer method, show that it always works for definite matrix pairs, and give its algorithm. At the end of the section we briefly describe the method of Zimmermann from [4] and [19], because it is closely related with the Falk–Langemeyer method. This relationship is also described.

The Falk–Langemeyer method solves problem (1) by constructing a sequence of “congruent” matrix pairs

(A^{(1)}, B^{(1)}), (A^{(2)}, B^{(2)}), …  (21)

where

A^{(1)} = A ,  B^{(1)} = B ,

A^{(k+1)} = F_{k}^{T} A^{(k)}F_{k} ,  B^{(k+1)} = F_{k}^{T} B^{(k)}F_{k} ,  k ≥ 1 . (22)

Note that the transformation (22) with nonsingular matrix F_{k} preserves the eigenvalues of the pair (A^{(k)}, B^{(k)}). This is a Jacobi–type method, hence the transformation matrices are chosen as nonsingular elementary plane matrices. An elementary plane matrix F = (f_{ij}) differs from the identity matrix only at the positions (l, l), (l, m), (m, l) and (m, m), where 1 ≤ l < m ≤ n. The matrix

F̂ = [ \begin{matrix}f_{ll} & f_{lm} \\ f_{ml} & f_{mm} \end{matrix} ]

is called (l, m)–restriction of the square matrix F = (f_{ij}).

For each k ≥ 1, the (l, m)–restriction of the matrix F_{k} has the form

F̂_{k} = [ \begin{matrix}1 & α_{k} \\ −β_{k} & 1 \end{matrix} ] ,  (23)

where real parameters α_{k} and β_{k} are chosen to satisfy the condition

a_{lm}^{(k+1)} = 0 ,  b_{lm}^{(k+1)} = 0 .  (24)

Here A^{(k)} = (a_{ij}^{(k)}) and B^{(k)} = (b_{ij}^{(k)}). Indices l and m are called pivot indices and the pair (l, m) is called pivot pair. As k varies the pivot pair also varies, hence l = l(k) and m = m(k). The transition from the pair (A^{(k)}, B^{(k)}) to the pair (A^{(k+1)}, B^{(k+1)}) is called the k–th step of the method. The manner in which we choose elements which are to be annihilated in the k–th step (or just the indices (l, m) of these elements) is called pivot strategy. The pivot strategy is cyclic if every sequence of N = n(n − 1)/2 successive pairs (l, m) contains all pairs (i, j), 1 ≤ i < j ≤ n. A sequence of N successive steps is referred to as a cycle. Two most common cyclic pivot strategies are the column–cyclic strategy and the row–cyclic strategy. The former is defined by the sequence of pairs

(1, 2), (1, 3), (2, 3), (1, 4), (2, 4), (3, 4), …, (1, n), (2, n), …, (n − 1, n),

and the latt er by the sequence of pairs

(1, 2), (1, 3), …, (1, n), (2, 3), (2, 4), …, (2, n), (3, 4), …, (n − 1, n).

These two strategies are also called serial strategies. Parallel cyclic strategies are cyclic strategies which enable simultaneous execution of approximately n/2 steps on parallel computers. These strategies have recently attracted considerable attention (see [12], [10]). We state the quadratic convergence results for serial and parallel strategies in Section 5.

Note that if the eigenvectors are needed, we must calculate the sequence of matrices F^{(1)}, F^{(2)}, …, where

F^{(1)} = I ,  F^{(k+1)} = F^{(k)}F_{k} ,  k ≥ 1 .  (25)

From the relations (22) and (25) we obtain for k ≥ 2

F^{(k)} = F_{1} ⋯ F_{k−1}
A^{(k)} = (F^{(k)})^{T} A^{(1)}F^{(k)} ,  B^{(k)} = (F^{(k)})^{T} B^{(1)}F^{(k)} .

We shall now derive one step of the algorithm. Note that only (l, m)– restrictions of the involved matrices are needed. Since (22) is the congruence transformation, the pairs (A^{(k)}, B^{(k)}) are definite for every k ≥ 1 and the pairs of the corresponding (l, m)–restrictions are definite as well.

Let (index k is omitted for the sake of simplicity)

F̂^{T} ÂF̂ = Â′ ,  F̂^{T} B̂F̂ = B̂′ ,

or respectively

[ \begin{matrix}1 & −β \\ α & 1 \end{matrix} ] [ \begin{matrix}a_{ll} & a_{lm} \\ a_{lm} & a_{mm} \end{matrix} ] [ \begin{matrix}1 & α \\ −β & 1 \end{matrix} ] = [ \begin{matrix}a′_{ll} & a′_{lm} \\ a′_{lm} & a′_{mm} \end{matrix} ] ,

[ \begin{matrix}1 & −β \\ α & 1 \end{matrix} ] [ \begin{matrix}b_{ll} & b_{lm} \\ b_{lm} & b_{mm} \end{matrix} ] [ \begin{matrix}1 & α \\ −β & 1 \end{matrix} ] = [ \begin{matrix}b′_{ll} & b′_{lm} \\ b′_{lm} & b′_{mm} \end{matrix} ] .  (26)

Condition (24) now reads a′_{lm} = b′_{lm} = 0. From the relation (26) the system in unknowns α and β is obtained

a′_{lm} = αa_{ll} + (1 − αβ)a_{lm} − βa_{mm} = 0 ,
b′_{lm} = αb_{ll} + (1 − αβ)b_{lm} − βb_{mm} = 0 .  (27)

Eliminating nonlinear terms in both equations we obtain

α = \frac{ℑ_{m}}{ν} ,  β = \frac{ℑ_{l}}{ν} ,  (28)

where ν is solution of the equation

ν^{2} − ℑ_{lm}ν − ℑ_{l}ℑ_{m} = 0  (29)

and

ℑ_{l} = a_{ll}b_{lm} − b_{ll}a_{lm} ,
ℑ_{m} = a_{mm}b_{lm} − b_{mm}a_{lm} ,
ℑ_{lm} = a_{ll}b_{mm} − a_{mm}b_{ll} .  (30)

Defining

ℑ = (ℑ_{lm})^{2} + 4ℑ_{l}ℑ_{m}  (31)

we obtain

ν_{±} = \frac{1}{2} sgn (ℑ_{lm})(| ℑ_{lm} | ±\sqrt{ℑ}) .

The algorithm is more stable if α and β are smaller in modulus, so we take

ν = ν_{+} = \frac{1}{2} sgn (ℑ_{lm})(| ℑ_{lm} | +\sqrt{ℑ}) .  (32)

From the above fromula we see that the necessary condition for carrying out this step is ℑ ≥ 0. Let us show that this condition is fulfilled in each step due to the definitness of the pairs ((A^{(k)},^{(k)} ), k ≥ 1). Proposition 3 Let the pair (A, B) be definite. Then the following holds:

(i) ℑ ≥ 0 , (ii) The following statments are equivalent:

(a) ℑ = 0 ,

(b) ℑ_{lm} = ℑ_{l} = ℑ_{m} = 0 ,

(c) There exist real constants s and t, |s| + |t| ≥ 0, such that

sÂ + tB̂ = O .

Proof: The proof can be found in [4] and [13] but for the completeness of exposition we present it below.

Using the relation (6) we can define the pair (A_{φ}, B_{φ}) such that the matrix B_{φ} is positive definite. Let us calculate the quantities (ℑ_{lm})_{φ}, (ℑ_{l})_{φ}, (ℑ_{m})_{φ} and (ℑ)_{φ} from the pair (A_{φ}, B_{φ}) using the relations (30) and (31). It is easy to verify that

ℑ_{l} = (ℑ_{l})_{φ},  ℑ_{m} = (ℑ_{m})_{φ},
ℑ_{lm} = (ℑ_{lm})_{φ},  ℑ = (ℑ)_{φ}.

Therefore, without loss of generality, we can assume that the matrix B from the pair (A, B) is positive definite. The statement (c) is now equivalent to the statement Â = cB̂, c ∈ R .

(i) With the notation

x = a_{ll}\sqrt{\frac{b_{mm}}{b_{ll}}} ,  y = a_{mm}\sqrt{\frac{b_{ll}}{b_{mm}}} ,  z = \frac{b_{lm}}{\sqrt{b_{ll}b_{mm}}} ,

the following identity holds:

ℑ = b_{ll}b_{mm}[(x − y)^{2} + 4(xz − a_{lm})(yz − a_{lm})].

Since the right side of the above relation is the square polynomial in a_{lm}, we have

ℑ = b_{ll}b_{mm}P_{2}(a_{lm}) ≥ b_{ll}b_{mm}P_{2} (\frac{x + y}{2} z)
= b_{ll}b_{mm}(x − y)^{2}(1 − z^{2}) = ℑ_{lm}^{2} (1 − \frac{b_{lm}^{2}}{b_{ll}b_{mm}} ) ≥ 0 .  (33)

In the last inequality we have used the assumption that the matrix B, and therefore the matrix B̂, is positive definite.

(ii) Let (a) hold. The relation (33) implies that ℑ_{lm} = 0. Matrix B is by the assumption positive definite. Therefore b_{ll} > 0, b_{mm} > 0, and the equality ℑ_{lm} = 0 can be written as

\frac{b_{mm}}{b_{ll}} a_{ll} = a_{mm}.

Using this relation we can write

ℑ_{m} = a_{mm}b_{lm} − b_{mm}a_{lm} = \frac{b_{mm}}{b_{ll}} (a_{ll}b_{lm} − b_{ll}a_{lm}) = \frac{b_{mm}}{b_{ll}} ℑ_{l} ,

or b_{ll}ℑ_{m} = b_{mm}ℑ_{l}. From the definition of ℑ, since ℑ_{lm} = 0, we conclude that ℑ_{l} = ℑ_{m} = 0. This gives (b).

Let (b) hold. Then

a_{mm} = b_{mm} \frac{a_{ll}}{b_{ll}} , a_{lm} = b_{lm} \frac{a_{ll}}{b_{ll}} .

Therefore, Â = cB̂, where

c = \frac{a_{ll}}{b_{ll}} = \frac{a_{mm}}{b_{mm}} ,

and (c) holds.

Let (c) hold. Then obviously (b) holds, and if (b) holds, then (a) holds.

Q.E.D.

Now we see that the Falk–Langemeyer method can be applied to all definite matrix pairs. Note that definitizing shifts are not used and need not to be known.

We have two special cases in the algorithm. If ℑ = 0, then the matrices Â and B̂ are proportional as shown in Proposition 3. Therefore, the two equations in the system (27) are linearly dependent and the system has a parametric solution in one of the following forms:

(α, β) = (\frac{c b_{mm} − b_{lm}}{b_{ll} − c b_{lm}} , c) ,  (α, β) = (\frac{c a_{mm} − a_{lm}}{a_{ll} − c a_{lm}} , c) ,
(α, β) = (c, \frac{c b_{ll} + b_{lm}}{b_{mm} + c b_{lm}} ) ,  (α, β) = (c, \frac{c a_{ll} + a_{lm}}{a_{mm} + c a_{lm}} ) ,

where c is real. For every c at least one of the quotients is well defined due to the definiteness of the pair (Â, B̂). It is best to set c = 0 to ensure that α_{k} and β_{k} tend to zero together with ε(A^{(k)}, B^{(k)}) as k → ∞ (see step (5a) in Algorithm 4). Setting c = 0 also reduces the operation count. This choice yields four possibilities for (α, β):

(−\frac{b_{lm}}{b_{ll}} , 0) ,  (−\frac{a_{lm}}{a_{ll}} , 0) ,  (34)
(0 , \frac{b_{lm}}{b_{mm}}) ,  (0 , \frac{a_{lm}}{a_{mm}}) .  (35)

Due to the definiteness of the pair (A, B), we have

|a_{ii}| + |b_{ii}| > 0,  i = 1, …, n ,  (36)

so at least one quotient is defined in each of the relations (34) and (35). In order to obtain better condition of the transformation matrix, we choose the relation in which the defined quotient has smaller absolute value. If both quotients in the chosen relation are defined, then they are equal, and for numerical reasons it is better to choose one in which the sum of squares of the numerator and the denominator is greater.

The second special case is when ℑ > 0 and ℑ_{lm} = 0 . This means that diagonals of the matrices Â and B̂ are proportional, while the matrices themselves are not. Then sgn (ℑ_{lm}) is not defined. Since ℑ_{l}ℑ_{m} > 0, we have sgn (ℑ_{l}) = sgn (ℑ_{m}) . Substituting sgn (ℑ_{lm}) with sgn (ℑ_{l}) in the equation (32) gives

ν = sgn (ℑ_{l})\sqrt{ℑ_{l}ℑ_{m}} .

Inserting this in the equation (28) gives, after simple calculation,

α = \sqrt{\frac{b_{mm}}{b_{ll}}} = \sqrt{\frac{a_{mm}}{a_{ll}}} ,  β = \frac{1}{α} .  (37)

The relation (36) implies that at least one of the quotients b_{mm}/b_{ll} and a_{mm}/a_{ll} is defined and different from zero. If both quotients are defined then they are equal and it is better to choose one in which the sum of squares of the numerator and the denominator is greater.

We can now define an algorithm of the method:

Algorithm 4 Definite matrix pair (A, B) is given.

(1) Set k = 1, A^{(1)} = A, B^{(1)} = B, F^{(1)} = I and choose the pivot strategy.

(2) Choose the pivot pair (l, m) = (l(k), m(k)) according to the strategy.

(3) If a_{lm}^{(k)} = b_{lm}^{(k)} = 0 , then set k = k + 1, A^{(k+1)} = A^{(k)}, B^{(k+1)} = B^{(k)}, F^{(k+1)} = F^{(k)} and go to step (2). Otherwise go to step (4).

(4) Calculate the quantities ℑ_{l}^{(k)}, ℑ_{m}^{(k)}, ℑ_{lm}^{(k)} and ℑ^{(k)} from formulas

ℑ_{l}^{(k)} = a_{ll}^{(k)}b_{lm}^{(k)} − b_{ll}^{(k)}a_{lm}^{(k)} ,  ℑ_{m}^{(k)} = a_{mm}^{(k)}b_{lm}^{(k)} − b_{mm}^{(k)}a_{lm}^{(k)} ,
ℑ_{lm}^{(k)} = a_{ll}^{(k)}b_{mm}^{(k)} − a_{mm}^{(k)}b_{ll}^{(k)} ,  ℑ^{(k)} = (ℑ_{lm}^{(k)})^{2} + 4ℑ_{l}^{(k)}ℑ_{m}^{(k)} .

(5)

(a) If ℑ^{(k)} = 0 perform the following steps: If | b_{ll}^{(k)} | ≥ | a_{ll}^{(k)} | , then set α_{k} = −b_{lm}^{(k)}/b_{ll}^{(k)} ;

otherwise set α_{k} = −a_{lm}^{(k)}/a_{ll}^{(k)} .

If | b_{mm}^{(k)} | ≥ | a_{mm}^{(k)} | , then set β_{k} = b_{lm}^{(k)}/b_{mm}^{(k)} ;

otherwise set β_{k} = a_{lm}^{(k)}/a_{mm}^{(k)} .

Finally, if | α_{k} | ≥ | β_{k} | , then set α_{k} = 0 ; otherwise set β_{k} = 0 .

(b) If ℑ^{(k)} > 0 perform the following steps:

(i) If ℑ_{lm}^{(k)} ≠ 0 , then calculate

ν_{k} = \frac{1}{2} sgn (ℑ_{lm}^{(k)})(| ℑ_{lm}^{(k)} | +\sqrt{ℑ^{(k)}}),

α_{k} = \frac{ℑ_{m}^{(k)}}{ν_{k}} , β_{k} = \frac{ℑ_{l}^{(k)}}{ν_{k}} .

(ii) If ℑ_{lm}^{(k)} = 0 , then, according to the relation (37), calculate

α_{k} = \sqrt{\frac{b_{mm}^{(k)}}{b_{ll}^{(k)}}} = \sqrt{\frac{a_{mm}^{(k)}}{a_{ll}^{(k)}}} , β_{k} = \frac{1}{α_{k}} .

If both quotients for α_{k} are defined, then choose one in which the sum of squares of the numarator and the denominator is greater.

(6) Perform the transformation

A^{(k+1)} = F_{k}^{T} A^{(k)}F_{k},  B^{(k+1)} = F_{k}^{T} B^{(k)}F_{k} ,  (38)
F^{(k+1)} = F^{(k)}F_{k} .  (39)

(7) Set k = k + 1 and move to step (2).

Since matrices A^{(k)}, B^{(k)}, A^{(k+1)} and B^{(k+1)} are symmetric, it is enough to store and to transform only upper triangles. In the transformation (38) only l−th and m−th row and column of the matrices A^{(k)} and B^{(k)} are changed and in the transformation (39) only l−th and m−th columns of the matrix F^{(k)} are changed. Note that the eigenvalues can be found without calculating the matrices F^{(k)} , k ≥ 1 , and therefore the trasformation (39) can be omitted. This reduces the operational count about fifty percent.

Stopping criteria of the infinite iterative procedure defined with this algorithm are described in Section 5.

From now on, the term “Falk–Langemeyer method” denotes the method described by Algorithm 4.

The Zimmerman method. We shall now relate the Falk–Langemeyer method with another method for solving the generalized eigenvalue problem. This method is due to K. Zimmermann who roughly described it in her thesis [19]. Later on, in his thesis [4], Hari derived its algorithm and proved its quadratic convergence.

The Zimmermann method is defined for symmetric matrix pairs (A, B) where matrix B is positive definite. We shall denote this fact as B > 0. At the beginning of the iterative procedure the initial pair (A, B) is normalized such that

A^{(1)} = DAD , B^{(1)} = DBD ,

where

D = diag (\frac{1}{\sqrt{b_{11}}} , …, \frac{1}{\sqrt{b_{nn}}}) .

Therefore, b_{ii}^{(1)} = 1, i = 1, …, n. The Zimmerman method constructs a sequence of pairs ((A^{(k)}, B^{(k)}), k ≥ 1) by the rule

A^{(k+1)} = Z_{k}^{T} A^{(k)}Z_{k} ,  B^{(k+1)} = Z_{k}^{T} B^{(k)}Z_{k} ,  k ≥ 1 .

The nonsingular matrices Z_{k} are chosen to preserve the units on the diagonal of B^{(k+1)} (automatic normalization at each step) and to annihilate the pivot elements. In [4] it is shown that for k ≥ 1 holds

Ẑ_{k} = \frac{1}{\sqrt{1 − (b_{lm}^{(k)})^{2}}} [ \begin{matrix}cos φ_{k} & sin φ_{k} \\ − sin ψ_{k} & cos ψ_{k} \end{matrix} ] ,

where

cos φ_{k} = cos θ_{k} + ξ_{k}(sin θ_{k} − η_{k} cos θ_{k}) ,
sin φ_{k} = sin θ_{k} − ξ_{k}(cos θ_{k} + η_{k} sin θ_{k}) ,
cos ψ_{k} = cos θ_{k} − ξ_{k}(sin θ_{k} + η_{k} cos θ_{k}) ,
sin ψ_{k} = sin θ_{k} + ξ_{k}(cos θ_{k} − η_{k} sin θ_{k}) ,

ξ_{k} = \frac{b_{lm}^{(k)}}{\sqrt{1 + b_{lm}^{(k)}} + \sqrt{1 − b_{lm}^{(k)}}} ,

η_{k} = \frac{b_{lm}^{(k)}}{(1 + \sqrt{1 + b_{lm}^{(k)}})(1 + \sqrt{1 − b_{lm}^{(k)}})} ,

tan 2θ_{k} = \frac{2a_{lm}^{(k)} − (a_{ll}^{(k)} + a_{mm}^{(k)})b_{lm}^{(k)}}{(a_{mm}^{(k)} − a_{ll}^{(k)})\sqrt{1 − (b_{lm}^{(k)})^{2}}} ,

−\frac{π}{4} ≤ θ_{k} ≤ \frac{π}{4} .

If a_{lm}^{(k)} = b_{lm}^{(k)} = 0 we set θ_{k} = 0. If the (l, m)–restrictions of A^{(k)} and B(k) are proportional and a_{lm}^{(k)} and b_{lm}^{(k)} are not both equal to zero, we set θ_{k} = \frac{π}{4}.

If the matrix B is not positive definite but the pair (A, B) is, then there exists a definitizing shift μ such that the matrix A − μB is positive definite. If this shift is known in advance, then the Zimmermann method can be applied to the pair (A, B) in the sense that each Z_{k} is computed from the pair (B^{(k)}, A^{(k)} − μB^{(k)}).

Although the Zimmermann method seems quite different from the Falk– Langemeyer method, the two methods are closely related. The following theorem gives precise formulation of this relationship. For this occasion only we assume that in step (5a) of Algorithm 4 (that is when ℑ^{(k)} = 0), parameters α_{k} and β_{k} are computed according to the formulae (37). For this version of the Falk–Langemeyer method holds:

Theorem 5 Let A and B be symmetric matrices of order n and let B be positive definite. Let the sequences ((A^{(k)}, B^{(k)}), k ≥ 1) and ((A^{(k)}′, B^{(k)}′), k ≥ 1) be generated from the starting pair (A, B) with the Falk–Langemeyer and the Zimmermann method, respectively. If the corresponding pivot strategies are the same, then

A^{(k)}′ = D^{(k)}A^{(k)}D^{(k)},  B^{(k)}′ = D^{(k)}B^{(k)}D^{(k)},  k ≥ 1 ,

where

D^{(k)} = diag (\frac{1}{\sqrt{b_{11}^{(k)}}} , …, \frac{1}{\sqrt{b_{nn}^{(k)}}} ) , k ≥ 1 .

Proof: The proof of this theorem is found in [4] Section 2.3. Q.E.D.

Let us suppose again that the matrix B is not positive definite while the pair (A, B) is, and that a positive definitizing shift μ is known in advance. Let us apply to the pair (A, B) the Zimmermann method in the sence mentioned above and the version of the Falk–Langemeyer method which we used in Theorem 5. It is easy to see that the parameters α_{k} and β_{k} from the Falk–Langemeyer method are invariant under the transformations (A, B) → (B, A − μB). Therefore, Theorem 5 holds in this case, as well, with

D^{(k)} = diag (\frac{1}{\sqrt{a_{11}^{(k)} − μ b_{11}^{(k)}}} , …, \frac{1}{\sqrt{a_{nn}^{(k)} − μ b_{nn}^{(k)}}} ) , k ≥ 1 .

We can conclude that if the starting pair is positive definite or the definitizing shift is known in advance, then the Falk–Langemmeyer (Zimmermann) method is the fast scaled (normalized) version of the Zimmermann (Falk– Langemmeyer) method.

4 Quadratic convergence

In this section we prove that the Falk–Langemeyer method is quadratically convergent if the starting definite pair has simple eigenvalues and the pivot strategy is cyclic. Definitizing shifts are not used and need not to be known. We first state the result about the quadratic convergence of the Zimmermann’s method, and show to what extent can this result be applied to the Falk–Langemeyer method if the matrix B is positive definite. Then we define the quadratic convergence for the Falk–Langemeyer method. In Subsection 4.1 we prove preliminary results which we use in the proof of the quadratic convergence of the Falk–Langemeyer method in Subsection 4.2.

The result about the quadratic convergence of Zimmermann method can be summarized as follows. Let the sequence ((A^{(k)}, B^{(k)}), k ≥ 1) be generated by the Zimmerman method from the pair (A, B), B > 0, and let ε_{k} = ε(A^{(k)}, B^{(k)}), where ε is defined with the relation (7). Note that ε_{k} is natural measure for convergence of the Zimmerman method since each matrix B^{(k)} has units along the diagonal.

We say that the Zimmerman method is quadratically convergent on the pair (A, B) if ε_{k} → 0 as k → ∞ and there exist a constant c_{0} > 0 and an integer r_{0} such that for r ≥ r_{0} holds

ε_{(r+1)N+1} ≤ c_{0} ε_{rN+1}^{2} .

Hence of special importance are conditions under which the above relation holds for r = 1. We call them asymptotic assumptions. Let

σ = spr (A, B) = max_{1≤i≤n} |λ_{i}| , δ = \frac{1}{3} min_{i≠j} | λ_{i} − λ_{j} | .

Theorem 6 Let the sequence ((A^{(k)}, B^{(k)}, k ≥ 1) be generated by the Zimmerman method from the starting pair (A, B), B > 0, and let the asymptotic assumptions

S(B^{(1)}) ≤ \frac{1}{2N} ,  2 \sqrt{1 + σ^{2}} ε_{1} < δ ,  (40)

hold. If the eigenvalues of the pair (A, B) are simple and the pivot strategy is cyclic, then

ε_{N+1} ≤ \sqrt{N(1 + σ^{2})} \frac{ε_{1}^{2}}{δ} .  (41)

Proof: The proof of this theorem is found in [4] Section 3.3. Q.E.D.

In Theorem 6 the term σ appears in the assumption (40) and in the assertion (41) because matrix B is not diagonal and matrix A is not normalized. From Theorem 5 we see that Theorem 6 holds for the Falk–Langemeyer method provided that the step (5a) of Algorithm 4 is appropriately changed, the matrix B is positive definite, and the pairs (A^{(k)}, B^{(k)}) generated by the Falk–Langemeyer method are normalized so that b_{ii}^{(k)} = 1, i = 1 … n, k ≥ 1.

In the rest of this section we prove that the Falk–Langemeyer method defined with Algorithm 4 is quadratically convergent on definite matrix pairs with simple eigenvalues if the pivot strategy is cyclic. We first have to define the measure for the quadratic convergence.

Let (A, B) be a definite pair. We shall use the measure ε̃ = ε̃(A, B) defined by

ε̃(A, B) = ε(Ã, B̃),

where Ã and B̃ are given by the relations (12),(5) and (4). The measure ε̃ enables us to use results of Corollary 2 and it takes into account the diagonal elements of matrices A and B. Note that the measure ε(A, B) is generally not the proper measure for almost diagonality of the pair (A, B) since it takes no account of the diagonals of matrices A and B.

Let the sequence of pairs

(A^{(1)}, B^{(1)}), (A^{(2)}, B^{(2)}), …  (42)

be generated by the Falk–Langemeyer method from the starting definite pair (A, B). For k ≥ 1 we set

ε̃_{k} = ε̃(A^{(k)}, B^{(k)}) = ε(Ã^{(k)}, B̃^{(k)}) ,  (43)
Ã^{(k)} = D_{k}A^{(k)}D_{k} ,  B̃^{(k)} = D_{k}B^{(k)}D_{k} ,  (44)
D_{k} = diag (\frac{1}{d_{1}^{(k)}} , …, \frac{1}{d_{n}^{(k)}} ) ,  (45)
d_{i}^{(k)} = ^{4} \sqrt{(a_{ii}^{(k)})^{2} + (b_{ii}^{(k)})^{2}} , i = 1, …, n .  (46)

From the relations (44), (45) and (46) we see that the pairs (Ã^{(k)}, B̃^{(k)}) are normalized in the sence that

(ã_{ii}^{(k)})^{2} + (b̃_{ii}^{(k)})^{2} = 1 ,  i = 1, …, n .  (47)

Definition 7 The Falk-Langemeyer method is quadratically convergent on the pair (A, B) if ε̃_{k} → 0 as k → ∞ and there exist a constant c_{0} > 0 and an integer r_{0} such that for r ≥ r_{0} holds

ε̃_{(r+1)N+1} ≤ c_{0}ε̃_{rN+1}^{2}.  (48)

From Definition 7 we see that ultimately ε̃_{k} decreases quadratically per cycle. At the end of Subsection 4.2 we shall show that the quadratic convergence implies the convergence of the sequence (42) towards the pair of diagonal matrices (D_{A}, D_{B}), where

D_{A} = diag (a_{1}, …, a_{n}),  D_{B} = diag (b_{1}, …, b_{n}) .  (49)

Here λ_{i} = [a_{i}, b_{i}], i = 1, …, n, are the eigenvalues of the pair (A, B). Finally, we shall show that ultimately the quadratic reduction of ε̃_{rN+1} implies the quadratic reduction of ε_{rN+1} and vice versa^{1}.

In order to be able to observe the measure ε̃ we must solve one more problem. The transformation matrices F_{k} are calculated from unnormalized pairs (A^{(k)}, B^{(k)}) and are therefore difficult to estimate. To solve this problem we shall observe the sequence obtained from the pair (A, B) with following process:

normalization, step of the method, normalization, step of the method,...

This sequence reads

( ˜\overline{A}^{(1)}, ˜ \overline{B}^{(1)}), (\overline{A}^{(2)}, \overline{B}^{(2)}), ( ˜\overline{A}^{(2)}, ˜ \overline{B}^{(2)}), (\overline{A}^{(3)}, \overline{B}^{(3)}), ( ˜\overline{A}^{(3)}, ˜ \overline{B}^{(3)}), …, (50)

where

( ˜\overline{A}^{(1)}, ˜ \overline{B}^{(1)}) = (Ã^{(1)}, B̃^{(1)}) ,  (51)

and for k ≥ 1 holds

\overline{A}^{(k+1)} = ˜ \overline{F} _{k}^{T} ˜ \overline{A}^{(k)} ˜\overline{F}_{k} ,  \overline{B}^{(k+1)} = ˜ \overline{F} _{k}^{T} ˜ \overline{B}^{(k)} ˜\overline{F}_{k} ,  (52)

1Here ε_{k} measures off–diagonal elements of the pairs from sequence (42) and should not be confused with the quantity used in connection with Zimmerman method.

˜\overline{A}^{(k+1)} = \overline{D}_{k+1}\overline{A}^{(k+1)}\overline{D}_{k+1} ,  ˜  \overline{B}^{(k+1)} = \overline{D}_{k+1}\overline{B}^{(k+1)}\overline{D}_{k+1} , (53)
\overline{D}_{k+1} = diag (\frac{1}{\overline{d}_{1}^{(k+1)}} , …, \frac{1}{\overline{d}_{n}^{(k+1)}} ) ,  (54)

\overline{d}_{i}^{(k+1)} = ^{4} \sqrt{(\overline{a}_{ii}^{(k+1)})^{2} + (\overline{b}_{ii}^{(k+1)})^{2}} ,  i = 1, …, n .  (55)

Of course, the sequences (42) and (50) are generated using the same pivot strategy. The matrices ˜ \overline{F}_{k} are calculated according to Algorithm 4, but

from the pairs ( ˜\overline{A}^{(k)}, ˜ \overline{B}^{(k)}). Since in the transition from ( ˜\overline{A}^{(k)}, ˜ \overline{B}^{(k)}) to (\overline{A}^{(k+1)}, \overline{B}^{(k+1)}) of all diagonal elements only those at positions (l, l) and (m, m) are being changed, we conclude that

\overline{d}_{i}^{(k+1)} = ^{4} \sqrt{(\overline{a}_{ii}^{(k+1)})^{2} + (\overline{b}_{ii}^{(k+1)})^{2}} = ^{4} \sqrt{(ã_{ii}^{(k)})^{2} + (b̃_{ii}^{(k)})^{2}} = 1, i = 1, …, n, i ≠ l, m .

(56)

We will now show that the operations of normalization and of carrying out one step of the algorithm commute. This is equivalent to showing that

˜\overline{A}^{(k)} = Ã^{(k)} and ˜ \overline{B}^{(k)} = B̃^{(k)} for k ≥ 1.

Let F̃_{k} be the transformation matrices obtained according to Algorithm 4 from the pairs (Ã^{(k)}, B̃^{(k)}), k ≥ 1. The following proposition shows that the matrices F_{k} and F̃_{k} are simply related.

Proposition 8 For k ≥ 1 holds F̃_{k} = D_{k}^{−1}F_{k}D_{k}.

Proof: Because of the relations (44), (45) and (46) we have

ˆÃ^{(k)} = [ \begin{matrix}\frac{a_{ll}^{(k)}}{(d_{l}^{(k)})^{2}} & \frac{a_{lm}^{(k)}}{d_{l}^{(k)}d_{m}^{(k)}} \\ \frac{a_{lm}^{(k)}}{d_{l}^{(k)}d_{m}^{(k)}} & \frac{a_{mm}^{(k)}}{(d_{m}^{(k)})^{2}} \end{matrix} ] ,  ˆ  B̃^{(k)} = [ \begin{matrix}\frac{b_{ll}^{(k)}}{(d_{l}^{(k)})^{2}} & \frac{b_{lm}^{(k)}}{d_{l}^{(k)}d_{m}^{(k)}} \\ \frac{b_{lm}^{(k)}}{d_{l}^{(k)}d_{m}^{(k)}} & \frac{b_{mm}^{(k)}}{(d_{m}^{(k)})^{2}} \end{matrix} ] . (57)

The assertion is now obtained by simply using the relation (57) in Algorithm 4 and calculating the matrix F̃_{k}. Q.E.D.

Proposition 9 For k ≥ 1 the following holds:

(i) ( ˜\overline{A}^{(k)}, ˜ \overline{B}^{(k)}) = (Ã^{(k)}, B̃^{(k)}) , (ii) D_{k} = D_{1}\overline{D}_{2}\overline{D}_{3} ⋯ \overline{D}_{k} .

Proof: The proof is by induction in respect to k.

(i) For k = 1 the assertion holds due to the relation (51). Suppose that the assertion holds for some k ≥ 1. This means that

˜\overline{A}^{(k)} = Ã^{(k)} ,  ˜  \overline{B}^{(k)} = B̃^{(k)} ,  ˜  \overline{F}_{k} = F̃_{k} .  (58)

From the relation (52) it follows that \overline{A}^{(k+1)} = ˜ \overline{F} _{k}^{T} ˜ \overline{A}^{(k)} ˜\overline{F}_{k} , which, because of the relation (58), implies that \overline{A}^{(k+1)} = F̃_{k}^{T} Ã^{(k)}F̃_{k}. Since the relation (44) and Proposition 8 imply

\overline{A}^{(k+1)} = D_{k}F_{k}^{T} D_{k}^{−1}D_{k}A^{(k)}D_{k}D_{k}^{−1}F_{k}D_{k} = D_{k}F_{k}^{T} A^{(k)}F_{k}D_{k}
= D_{k}A^{(k+1)}D_{k} ,  (59)

we conclude that normalizations of the matrices A^{(k+1)} and \overline{A}^{(k+1)} give the

same matrix. Now we use the same argument to show that ˜ \overline{B}^{(k+1)} = B̃^{(k+1)} for k ≥ 1 and to prove (i).

(ii) For k = 1 the assertion is trivially fulfilled. Let the assertion hold for some k ≥ 1 . From the relations (59), (53) and the assertion (i) we obtain

\overline{D}_{k+1}D_{k}A^{(k+1)}D_{k}\overline{D}_{k+1} = \overline{D}_{k+1}\overline{A}^{(k+1)}\overline{D}_{k+1} =
= ˜ \overline{A}^{(k+1)} = Ã^{(k+1)} = D_{k+1}A^{(k+1)}D_{k+1} .

It is obvious that D_{k+1} = D_{k}\overline{D}_{k+1} and inserting the induction assumption we conclude that (ii) holds. Q.E.D.

From Proposition 9 we see that the relations (50), (52) and (53) can be written as

(Ã^{(1)}, B̃^{(1)}), (\overline{A}^{(2)}, \overline{B}^{(2)}), (Ã^{(2)}, B̃^{(2)}), (\overline{A}^{(3)}, \overline{B}^{(3)}), (Ã^{(3)}, B̃^{(3)}), … (60)
\overline{A}^{(k+1)} = F̃_{k}^{T} Ã^{(k)}F̃_{k} ,  \overline{B}^{(k+1)} = F̃_{k}^{T} B̃^{(k)}F̃_{k} ,  (61)
Ã^{(k+1)} = \overline{D}_{k+1}\overline{A}^{(k+1)}\overline{D}_{k+1} ,  B̃^{(k+1)} = \overline{D}_{k+1}\overline{B}^{(k+1)}\overline{D}_{k+1} . (62)

The relations (60), (61), (62), (54) and (55) define the normalized Falk– Langemeyer method. We use the normalized method only as an aid to obtain information about the quantity ε̃_{k}. æ

4.1 Preliminaries

Here we define asymptotic assumptions and prove several lemmas which are used later in the proof of the quadratic convergence of the Falk–Langemeyer method. The quadratic convergence proof is based on the idea of Wilkinson (see [18]) which consists in estimating the growth of already annihilated elements in the current cycle. To this end we must estimate the transformation parametars α̃_{k} and β̃_{k} and also the growth of all off-diagonal elements in the current cycle. These two tasks are solved in Lemma 11, Lemma 13, Lemma 14 and Lemma 15. Lemma 10 gives us two numeric relations which are used in the proof. Lemma 11 and Lemma 12 estimate the transformation parametars α̃_{k} and β̃_{k}, and the measure ε̃_{k} in one step. Lemma 13, Lemma 14 and Lemma 15 estimate the growth of α̃_{k}, β̃_{k} and ε̃_{k} during N consecutive steps. Lemma 15 is the most important for the proof of the quadratic convergence. In this subsection we do not assume that the pivot strategy is cyclic. Therefore the results of this subsection hold for any pivot strategy. However, if the pivot strategy is cyclic, then Lemma 13, Lemma 14 and Lemma 15 explain the behaviour of α̃_{k}, β̃_{k} and ε̃_{k} during one cycle.

As we said in Section 1, the quadratic convergence can always be expected if the eigenvalues of problem (1) are simple. We will therefore use two quadratic convergence assumptions:

(A1) The eigenvalues of the pair (A, B) are simple, i. e.

p = n ≥ 3 .

(A2) The pair (A, B) is almost diagonal, i. e.

\frac{ε̃_{1}}{δ} < \frac{1}{2N} .

Asymptotic assumption (A2) is similar to the assumptions used in Theorem 6 and in convergence results of some other Jacobi–type methods (see [4], [1]). Assumption (A1) implies

N ≥ 3  (63)

and

ε_{k} = τ_{k} ,  ε̃_{k} = τ̃_{k} ,  k ≥ 1 ,  (64)

where τ_{k} = τ(A^{(k)}, B^{(k)}) and τ̃_{k} = τ(Ã^{(k)}, B̃^{(k)}). We shall use the notation

ã_{k} = |ã_{lm}^{(k)}| ,  b̃_{k} = |b̃_{lm}^{(k)}| ,  k ≥ 1 .  (65)

Lemma 10 Let r be an integer such that r ≥ 3 and let x be a nonnegative real number satisfiying 2xr < 1. Then the following inequalities hold:

(1 − x)^{−r} ≤ 1 + \frac{12}{7} · r · x , (1 + x)^{r} ≤ 1 + \frac{4}{3} · r · x .

Proof: The proof of this lemma is elementary and can be found in [4]. Q.E.D.

The following lemma shows how are the transformation parameters α̃_{k} and β̃_{k} from matrices F̃_{k} bounded with ϵ̃_{k}.

Lemma 11 Let the assumption (A1) hold. If for some k ≥ 1 holds

ε̃_{k} < \frac{2}{3N} δ ,  (66)

then

max \{|α̃_{k}|, |β̃_{k}|\} ≤ 0.34 · \frac{\sqrt{(ã_{k})^{2} + (b̃_{k})^{2}}}{δ} .  (67)

Proof: Suppose that for some k ≥ 1 the relation (66) holds. Then Theorem 1 and Corollary 2 hold for the pair (Ã^{(k)}, B̃^{(k)}) as well. The assumption (A1) and the relations (63), (64) and (18) imply that there exists an ordering of the eigenvalues of the pair (A, B)^{2} such that

χ(λ_{i}, [ã_{ii}^{(k)}, b̃_{ii}^{(k)}]) ≤ \frac{ε̃_{k}^{2}}{2δ} < \frac{4δ^{2}}{9N^{2}} · \frac{1}{2δ} < \frac{2}{81} · δ < 0.025 · δ ,
i = 1, …, n . (68)

Applying twice the triangle inequality and using the definition (11) and the relation (68), we obtain

| ℑ̃_{lm}^{(k)} | = | ã_{ll}^{(k)}b̃_{mm}^{k)} − ã_{mm}^{(k)}b̃_{ll}^{(k)} | = χ([ã_{ll}^{(k)}, b̃_{ll}^{(k)}], [ã_{mm}^{(k)}, b̃_{mm}^{(k)}])
≥ χ(λ_{l}, λ_{m}) − χ(λ_{l}, [ã_{ll}^{(k)}, b̃_{ll}^{(k)}]) − χ(λ_{m}, [ã_{mm}^{(k)}, b̃_{mm}^{(k)}])
> 3 · δ − 2 · 0.025 · δ = 2.95 · δ .  (69)

2Since p = n, the eigenvalues can be ordered so that the matrix P from Theorem 1 and Corollary 2 is identity matrix.

It is obvious that |ℑ̃_{lm}^{(k)}| ≠ 0. This excludes cases (5a) and (5bii) of Algorithm 3. Therefore, we have

max\{|α̃_{k}|, |β̃_{k}|\} ≤ \frac{2}{|ℑ̃_{lm}^{(k)}| + \sqrt{ℑ̃^{(k)}}} · max\{|ℑ̃_{l}^{(k)}| , |ℑ̃_{m}^{(k)}|\} .  (70)

From the Cauchy–Schwarz inequality and the relations (47) and (65) we have

|ℑ̃_{l}^{(k)}| = |ã_{ll}^{(k)}b̃_{lm}^{(k)} − b̃_{ll}^{(k)}ã_{lm}^{(k)}| ≤ \sqrt{(ã_{ll}^{(k)})^{2} + (b̃_{ll}^{(k)})^{2}}\sqrt{(ã_{lm}^{(k)})^{2} + (b̃_{lm}^{(k)})^{2}}
= \sqrt{(ã_{k})^{2} + (b̃_{k})^{2}} .

The same estimate holds for ℑ̃_{m}^{(k)} and therefore

max\{|ℑ̃_{l}^{(k)}| , |ℑ̃_{m}^{(k)}|\} ≤ \sqrt{(ã_{k})^{2} + (b̃_{k})^{2}} .  (71)

Since

(ã_{k})^{2} + (b̃_{k})^{2} ≤ \frac{1}{2} · ε̃_{k}^{2} ,  (72)

the relations (69), (71),and (72) imply

\sqrt{ℑ̃^{(k)}} = \sqrt{(ℑ̃_{lm}^{(k)})^{2} + 4 ℑ̃_{l}^{(k)}ℑ̃_{m}^{(k)}} ≥ \sqrt{(2.95 δ)^{2} − 4 · \frac{ε̃_{k}^{2}}{2}}
≥ \sqrt{(2.95 δ)^{2} − 2 · \frac{4}{9N^{2}} δ^{2}} ≥ 2.933 δ .

The assertion (67) now follows from the relations (70), (69), (71) and the above relation. Q.E.D.

The following Lemma gives the relation between ε̃_{k} and ε̃_{k+1}. It is used later in the proof of of Lemma 15.

Lemma 12 Let the assumption (A1) hold. If for some k ≥ 1 the relation (66) holds, then

ε̃_{k+1}^{2} ≤ \frac{1 + 0.494 · \frac{_{ε_{k}}}{δ}}{1 − 0.077 · \frac{_{ε_{k}}}{δ}} [ε̃_{k}^{2} − 2 (ã_{k}^{2} + b̃_{k}^{2})] ,  (73)

Proof: Suppose that the relation (66) holds for some k ≥ 1. The relation (62), together with the definition of ε̃_{k}, implies

ε̃_{k+1}^{2} = S^{2}(\overline{D}_{k+1}\overline{A}^{(k+1)}\overline{D}_{k+1}) + S^{2}(\overline{D}_{k+1}\overline{B}^{(k+1)}\overline{D}_{k+1}) .  (74)

If

m_{k+1} = min\{\overline{d}_{1}^{(k+1)}, …, \overline{d}_{n}^{(k+1)}\} ,

then the relation (74) implies

ε̃_{k+1}^{2} ≤ S^{2}(\frac{1}{(m_{k+1})^{2}} \overline{A}^{(k+1)}) + S^{2}(\frac{1}{(m_{k+1})^{2}} \overline{B}^{(k+1)}) = \frac{1}{(m_{k+1})^{4}} · \overline{ε}_{k+1}^{2} . (75)

Let us define vectors

ã^{l} = (ã_{l,1}^{(k)}, ã_{l,2}^{(k)}, …, ã_{l,l−1}^{(k)}, ã_{l,l+1}^{(k)}, …, ã_{l,m−1}^{(k)}, ã_{l,m+1}^{(k)}, …, ã_{l,n}^{(k)}) ,
ã^{m} = (ã_{m,1}^{(k)}, ã_{m,2}^{(k)}, …, ã_{m,l−1}^{(k)}, ã_{m,l+1}^{(k)}, …, ã_{m,m−1}^{(k)}, ã_{m,m+1}^{(k)}, …, ã_{m,n}^{(k)}) ,
ã_{l}^{T} = (ã_{1,l}^{(k)}, ã_{2,l}^{(k)}, …, ã_{l−1,l}^{(k)}, ã_{l+1,l}^{(k)}, …, ã_{m−1,l}^{(k)}, ã_{m+1,l}^{(k)}, …, ã_{n,l}^{(k)}) ,
ã_{m}^{T} = (ã_{1,m}^{(k)}, ã_{2,m}^{(k)}, …, ã_{l−1,m}^{(k)}, ã_{l+1,m}^{(k)}, …, ã_{m−1,m}^{(k)}, ã_{m+1,m}^{(k)}, …, ã_{n,m}^{(k)}) ,

where generally a^{T} denotes the transposed vector a. Let \overline{a}^{l}, \overline{a}^{m}, \overline{a}_{l} and \overline{a}_{m} be row and column vectors defined in the same way, but from elements of the matrix \overline{A}^{(k+1)}. Relation (61) implies that

[ \begin{matrix}\overline{a} & l \\ \overline{a} & m \end{matrix} ] = ˆ F̃_{k}^{T} [ ã^{l} ; ã^{m} ] , [\overline{a}_{l}, \overline{a}_{m}] = [ã_{l}, ã_{m}] ˆF̃_{k} .

Therefore,

∥[ \begin{matrix}\overline{a} & l \\ \overline{a} & m \end{matrix} ]∥^{2} ≤ ∥ ˆ F̃_{k}^{T} ∥_{2}^{2} ∥[ ã^{l} ; ã^{m} ]∥^{2} , ∥ [\overline{a}_{l}, \overline{a}_{m}] ∥^{2} ≤ ∥ ˆ F̃_{k} ∥_{2}^{2} ∥ [ã_{l}, ã_{m}] ∥^{2} .

The off–diagonal elements of the matrix Ã^{(k)} which are changed in the transformation (61), are exactly the elements of vectors ã^{l}, ã^{m}, ã_{l} and ã_{m} with the

exception of ã_{lm}^{(k)} and ã_{ml}^{(k)} which are annihilated. Since ∥ ˆ F̃_{k}^{T} ∥_{2} = ∥ ˆ F̃_{k} ∥_{2}, we conclude that

S^{2}(\overline{A}^{(k+1)}) ≤ S^{2}(Ã^{(k)})+(∥ ˆ F̃_{k} ∥_{2}^{2} −1)(∥ ã^{l} ∥^{2} + ∥ ã^{m} ∥^{2} + ∥ ã_{l} ∥^{2} + ∥ ã_{m} ∥^{2})−2ã_{k}^{2} .

Since ∥ ˆ F̃_{k} ∥_{2}≥ 1 (see further in the proof), we conclude that

S^{2}(\overline{A}^{(k+1)}) ≤ ∥ ˆ F̃_{k} ∥_{2}^{2} (S^{2}(Ã^{(k)}) − 2ã_{k}^{2}) .

By applying the similar analysis to matrix B̃^{(k)}, we obtain

S^{2}(\overline{B}^{(k+1)}) ≤ ∥ ˆ F̃_{k} ∥_{2}^{2} (S^{2}(B̃^{(k)}) − 2b̃_{k}^{2}) .

Adding two previous inequalities and using the definitions of \overline{ε}_{k+1} and ε̃_{k}, gives

\overline{ε}_{k+1}^{2} ≤ ∥ ˆ F̃_{k} ∥_{2}^{2} [ε̃_{k}^{2} − 2(ã_{k}^{2} + b̃_{k}^{2})] .

Inserting this inequality into relation (75), we obtain

ε̃_{k+1}^{2} ≤ \frac{∥ ˆ F̃_{k} ∥_{2}^{2}}{(m_{k+1})^{4}} · [ε̃_{k}^{2} − 2(ã_{k}^{2} + b̃_{k}^{2})] .  (76)

To complete the proof we must find the upper bound for ∥ ˆ F̃_{k} ∥_{2}^{2} and the lower bound for m_{k+1}.

The relation (56) implies that

m_{k+1} = min\{1, \overline{d}_{l}^{(k+1)}, \overline{d}_{m}^{(k+1)}\} .  (77)

Relation (61) implies that

\overline{a}_{ll}^{(k+1)} = ã_{ll}^{(k)} − 2β̃_{k}ã_{lm}^{(k)} + β̃_{k}^{2} ã_{mm}^{(k)} ,  \overline{b}_{ll}^{(k+1)} = b̃_{ll}^{(k)} − 2β̃_{k}b̃_{lm}^{(k)} + β̃_{k}^{2} b̃_{mm}^{(k)} ,
\overline{a}_{mm}^{(k+1)} = α̃_{k}^{2} ã_{ll}^{(k)} + 2α̃_{k}ã_{lm}^{(k)} + ã_{mm}^{(k)} ,  \overline{b}_{mm}^{(k+1)} = α̃_{k}^{2} b̃_{ll}^{(k)} + 2α̃_{k}b̃_{lm}^{(k)} + b̃_{mm}^{(k)} .

Therefore

(\overline{d}_{l}^{(k+1)})^{4} = (\overline{a}_{ll}^{(k+1)})^{2} + (\overline{b}_{ll}^{(k+1)})^{2} = 1 + 4 β_{k}^{2} [(ã_{lm}^{(k)})^{2} + (b̃_{lm}^{(k)})^{2}] +
+ β̃_{k}^{4} − 4 β̃_{k}(ã_{ll}^{(k)}ã_{lm}^{(k)} + b̃_{ll}^{(k)}b̃_{lm}^{(k)}) − 4 β̃_{k}^{3}(ã_{mm}^{(k)}ã_{lm}^{(k)} + b̃_{mm}^{(k)}b̃_{lm}^{(k)}) +
+ 2β̃_{k}^{2}(ã_{ll}^{(k)}ã_{mm}^{(k)} + b̃_{ll}^{(k)}b̃_{mm}^{(k)})
≥ 1 − 4 | β̃_{k} || ã_{ll}^{(k)}ã_{lm}^{(k)} + b̃_{ll}^{(k)}b̃_{lm}^{(k)} | −4 | β̃_{k} |^{3} | ã_{mm}^{(k)}ã_{lm}^{(k)} + b̃_{mm}^{(k)}b̃_{lm}^{(k)} | −
− 2 β̃_{k}^{2} | ã_{ll}^{(k)}ã_{mm}^{(k)} + b̃_{ll}^{(k)}b̃_{mm}^{(k)} | .  (78)

Using the relation (47) and the Cauchy–Schwarz inequality in the relation (78), we obtain

(\overline{d}_{l}^{(k+1)})^{4} ≥ 1 − 4 | β̃_{k} | (1+ | β̃_{k} |^{2})\sqrt{(ã_{lm}^{(k)})^{2} + (b̃_{lm}^{(k)})^{2}} − 2β̃_{k}^{2} .  (79)

Using similar argument we obtain

(\overline{d}_{m}^{(k+1)})^{4} ≥ 1 − 4 | α̃_{k} | (1+ | α̃_{k} |^{2})\sqrt{(ã_{lm}^{(k)})^{2} + (b̃_{lm}^{(k)})^{2}} − 2α̃_{k}^{2} .  (80)

Relations (77), (79), (80), (72) and Lemma 11 now imply

m_{k+1}^{4} ≥ 1 − 4 · \frac{0.34}{\sqrt{2}} · \frac{ε̃_{k}}{δ} · [ 1 + (\frac{0.34}{\sqrt{2}} · \frac{ε̃_{k}}{δ} )^{2}] · \frac{ε̃_{k}}{\sqrt{2}} − 2 · (\frac{0.34}{\sqrt{2}} · \frac{ε̃_{k}}{δ} )^{2} . (81)

Since χ is chordal metric, from the relation (11) we see that δ ≤ 1/3. The relations (66) and (63) therefore imply

ε̃_{k} < \frac{2}{9N} < \frac{2}{27} .  (82)

Inserting the relation (82) and the assumption (66) into the relation (81) we obtain

m_{k+1}^{4} > 1 − \frac{0.34}{\sqrt{2}} · \frac{ε̃_{k}}{δ} [ 4 · (1 + (\frac{0.34}{\sqrt{2}} · \frac{2}{3N} )^{2}) · \frac{2}{27 \sqrt{2}} + 0.34 · \sqrt{2} · \frac{2}{3N} ] .

Finally, taking into account that N ≥ 3 we obtain

m_{k+1}^{4} ≥ 1 − 0.077 · \frac{ε̃_{k}}{δ}  (83)

We shall now estimate ∥ ˆ F̃_{k} ∥_{2}^{2}. Since

∥ ˆ F̃_{k} ∥_{2}^{2} ≤ ∥ ˆ F̃_{k} ∥_{1} · ∥ ˆ F̃_{k} ∥_{∞} ,

where ∥A∥_{1} = max_{j} ∑_{i} |a_{ij}|, ∥A∥_{∞} = max_{i} ∑_{j} |a_{ij}| for A = (a_{ij}), we obtain

∥ ˆ F̃_{k} ∥_{2}^{2} ≤ (1 + max{|α̃_{k}|, |β̃_{k}|})^{2} .

Lemma 11 and the relation (66) now imply

∥ ˆ F̃_{k} ∥_{2}^{2} ≤ (1 + \frac{0.34}{\sqrt{2}} \frac{ε̃_{k}}{δ} )^{2} ≤ 1 + \frac{0.34}{\sqrt{2}} \frac{ε̃_{k}}{δ} (2 + \frac{0.34}{\sqrt{2}} \frac{2}{3N} )
≤ 1 + 0.494 \frac{ε̃_{k}}{δ} .  (84)

The relation (73) now follows from the relations (76), (83) and (84). Q.E.D.æ

We shall now prove that if the assumptions (A1) and (A2) hold, then Lemma 11 and Lemma 12 hold during N consecutive steps.

Lemma 13 Let the asymptotic assumptions (A1) and (A2) hold. Then for each k ∈ {1, …, N} holds

ε̃_{k} ≤ \frac{1}{1 − 0.3 · (k − 1)\frac{_{ε_{1}}}{δ}} · ε̃_{1} , \frac{ε̃_{k}}{δ} < \frac{2}{3N} .

Proof: The proof is by induction. For k = 1 lemma is trivially fulfilled. Suppose that lemma holds for some k ∈ {1, …, N − 1}. From the second inequality in the induction assumption we conclude that, for the chosen k, Lemma 11 and Lemma 12 hold. From Lemma 12 it follows that

ε̃_{k+1}^{2} ≤ \frac{1 + 0.494 · \frac{_{ε_{k}}}{δ}}{1 − 0.077 · \frac{_{ε_{k}}}{δ}} · ε̃_{k}^{2} ≤ \frac{1}{(1 − 0.494\frac{_{ε_{k}}}{δ})(1 − 0.077\frac{_{ε_{k}}}{δ})}ε̃_{k}^{2}
≤ \frac{1}{(1 − 0.3\frac{_{ε_{k}}}{δ})^{2}} ε̃_{k}^{2} .  (85)

Hence

ε̃_{k+1} ≤ \frac{1}{1 − 0.3 · \frac{_{ε_{k}}}{δ}} · ε̃_{k} .

Inserting the induction assumption in this inequality we obtain

ε̃_{k+1} ≤ \frac{1}{1 − 0.3_{\frac{1}{1−0.3(k−1)\frac{_{ε_{1}}}{δ}}} \frac{_{ε_{1}}}{δ}} · \frac{1}{1 − 0.3(k − 1)\frac{_{ε_{1}}}{δ}} · ε̃_{1}

≤ \frac{1}{1 − 0.3(k − 1)\frac{_{ε_{1}}}{δ} − 0.3\frac{_{ε_{1}}}{δ}} · ε̃_{1} = \frac{1}{1 − 0.3 · k · \frac{_{ε_{1}}}{δ}} · ε̃_{1} ,

and the first assertion of the lemma is proved. From this assertion for k + 1, because of the asymptotic assumption (A2), we now have

\frac{ε̃_{k+1}}{δ} ≤ \frac{1}{1 − 0.3 · k · \frac{_{ε_{1}}}{δ}} · \frac{ε̃_{1}}{δ} < \frac{1}{1 − 0.3 (N − 1)\frac{1}{2N}} · \frac{1}{2N} = \frac{1}{2 − 0.3} · \frac{1}{N} < \frac{2}{3N}

which completes the proof. Q.E.D.

Lemma 14 If the asymptotic assumptions (A1) and (A2) hold, then the assertions (67) and (73) of Lemma 11 and Lemma 12 hold for every k ∈ {1, …, N}.

Proof: The assertion follows imidiately from second assertion of Lemma 13.

Q.E.D.

The next lemma explains behaviour of S(Ã^{(k)}), S(B̃^{(k)}) and ε̃_{k}, and of the transformation parameters α̃_{k} and β̃_{k} during N consecutive steps. Let us define the quantity

c_{N} = \frac{1 + 0.494 · \frac{2}{3N}}{1 − 0.077 · \frac{2}{3N}} .  (86)

Lemma 15 Let the asymptotic assumptions (A1) and (A2) hold. Then:

(i) For k = 1, …, N holds

[ S^{2}(Ã^{(k+1)}) ; S^{2}(B̃^{(k+1)}) ; ε̃_{k+1}^{2} ] ≤ (c_{N})^{k} [ S^{2}(Ã^{(1)}) ; S^{2}(B̃^{(1)}) ; ε̃_{1}^{2} ] ≤ 1.566 [ S^{2}(Ã^{(1)}) ; S^{2}(B̃^{(1)}) ; ε̃_{1}^{2} ] .

(ii) For any choice ω̃_{k} ∈ {α̃_{k}, β̃_{k}}, 1 ≤ k ≤ N, holds

∑_{k=1}^{N} ω̃_{k}^{2} ≤ 0.426 · \frac{ε̃_{1}^{2}}{δ^{2}} .

Proof: (i) Because of Lemma 12 and Lemma 14, for k = 1, …, N holds

ε̃_{k+1}^{2} ≤ c_{N}(ε̃_{k}^{2} − 2 (ã_{k}^{2} + b̃_{k}^{2}))
≤ c_{N}{c_{N}[ ε̃_{k−1}^{2} − 2 (ã_{k−1}^{2} + b̃_{k−1}^{2})] − 2 (ã_{k}^{2} + b̃_{k}^{2})}
≤ … ≤ (c_{N})^{k}ε̃_{1}^{2} − 2 ∑_{j=1}^{k}(c_{N})^{k−j+1}(ã_{j}^{2} + b̃_{j}^{2}) .  (87)

From the relation (87) immediately follows

ε̃_{k+1}^{2} ≤ (c_{N})^{k}ε̃_{1}^{2} ≤ (c_{N})^{N} ε̃_{1}^{2} ,  k = 1, …, N .  (88)

Using Lemma 10 we obtain

(c_{N})^{N} < (1 + \frac{4}{3} · 0.494 · \frac{2}{3N} · N)(1 + \frac{12}{7} · 0.077 · \frac{2}{3N} · N) < 1.566 . (89)

Inserting this inequality into relation (88), we obtain

ε̃_{k+1}^{2} ≤ 1.566 · ε̃_{1}^{2} , k = 1, …, N .

From the proof of Lemma 12 we se that the above estimates hold for the quantities S^{2}(Ã^{(k+1)}) and S^{2}(B̃^{(k+1)}), as well. Therefore (i) is proved.

(ii) Since c_{N} > 1 , from the relation (87) for k = N we have

ε̃_{N+1}^{2} ≤ (c_{N})^{N} ε̃_{1}^{2} − 2 ∑_{k=1}^{N}(ã_{k}^{2} + b̃_{k}^{2}) .

Since ε̃_{N+1}^{2} ≥ 0, this inequality implies

∑_{k=1}^{N}(ã_{k}^{2} + b̃_{k}^{2}) ≤ \frac{1}{2} · (c_{N})^{N} ε̃_{1}^{2} ≤ 0.783 · ε̃_{1}^{2} .

The above inequality together with Lemma 11 and Lemma 14 imply

∑_{k=1}^{N} ω̃_{k}^{2} ≤ ∑_{k=1}^{N} max\{α̃_{k}^{2}, β̃_{k}^{2}\} ≤ ∑_{k=1}^{N} 0.34^{2} · (ã_{k}^{2} + b̃_{k}^{2}) · \frac{1}{δ^{2}}
≤ 0.1156 · 0.783 · \frac{ε̃_{1}^{2}}{δ^{2}} ≤ 0.091 · \frac{ε̃_{1}^{2}}{δ^{2}}

and the lemma is proved. Q.E.D.

4.2 The proof

Here we prove that the Falk–Langemeyer method is quadratically convergent if the assumptions (A1) and (A2) are fulfilled and the pivot strategy is cyclic. Then we prove that the quadratic convergence implies the convergence of the sequence of pairs (42) towards the pair of diagonal matrices. At the end we prove that the measures ε̃_{k} and ε_{k} are equivalent in the sense that ultimately the quadratic reduction of ε̃_{kN+1} implies the quadratic reduction of ε_{kN+1} and vice versa.

We can now prove our paper’s central theorem.

Theorem 16 Let the asymptotic assumptions (A1) and (A2) hold and let the sequence ((A^{(k)}, B^{(k)}), k ≥ 1) be generated with the Falk–Langemeyer method from the pair (A, B). Then for any cyclic strategy holds

ε̃_{N+1} ≤ \sqrt{N} · \frac{ε̃_{1}^{2}}{δ} .

Proof: Let us fix some k ∈ {1, …, N}. Then the pivot pair (l, m) is also fixed. We want to know what happens with the element on this position till the end of cycle. Therefore, we will observe the elements ã_{lm}^{(r)}, r = k + 1, …, N. We know that ã_{lm}^{(k+1)} = 0 and that the elements ã_{lm}^{(r)} actually change at most 2(n − 2) times. Let r_{1}, …, r_{s}, s ≤ 2n − 4, denote those values of r for which ã_{lm}^{(r)} changes in the r–th step. Let us introduce the notation:

\overline{z}_{i} = \overline{a}_{lm}^{(r_{i}+1)},  z̃_{i} = ã_{lm}^{(r_{i}+1)} ,
\overline{d}_{j}^{(i)} = ^{4} \sqrt{(\overline{a}_{jj}^{(r_{i}+1)})^{2} + \overline{b}_{jj}^{(r_{i}+1)})^{2}} ,  j ∈ \{l, m\} ,
\overline{m}_{lm}^{(i)} = min\{\overline{d}_{l}^{(i)}, \overline{d}_{m}^{(i)}\} ,  d_{N} = \sqrt{1 − 0.077 · \frac{2}{3N}} .  (90)

Performing the r_{1}–th step according to Algorithm 4, gives

\overline{z}_{1} = (0 · 1 ± ã^{(r_{1})}ω̃_{r_{1}}) ,

where ω̃_{r_{1}} ∈ {α̃_{r_{1}}, β̃_{r_{1}}} and ã^{(r_{1})} is some off–diagonal element of the matrix Ã^{(r_{1})}. Since

z̃_{i} = \begin{matrix}\overline{z} & i \\ \overline{d} & {(i) \atop l} & \overline{d} & {(i) \atop m} \end{matrix} ,  i = 1, …, s ,  (91)

from the relations (90), (83) and Lemma 12 follows that

| z̃_{1} | ≤ \frac{1}{(\overline{m}_{lm}^{(i)})^{2}} | ã^{(r_{1})} || ω̃_{r_{1}} | ≤ \frac{1}{d_{N}} | ã^{(r_{1})} || ω̃_{r_{1}} | .  (92)

Further, in the r_{2}–th step, we have

\overline{z}_{2} = (1 · z̃_{1} ± ã^{(r_{2})}ω̃_{r_{2}}) ,  (93)

where ω̃_{r_{2}} ∈ {α̃_{r_{2}}, β̃_{r_{2}}}, and ã^{(r_{2})} is some off–diagonal element of the matrix Ã^{(r_{2})}. The relations (93), (92) and (91) imply

| z̃_{2} | ≤ \frac{1}{d_{N}} · (\frac{1}{d_{N}} | ã^{(r_{1})} || ω̃_{r_{1}} | + | ã^{(r_{2})} || ω̃_{r_{2}} |) .

By induction we obtain

| z̃_{j} | ≤ ∑_{i=1}^{j} \frac{1}{(d_{N})^{j−i+1}} | ã^{(r_{i})} || ω̃_{r_{i}} | ,  j = 1, …, s .  (94)

For k = 1, …, N + 1 following notation is introduced:

Ã^{(k)} = D̃_{A}^{(k)} + Ẽ^{(k)} ,  D̃_{A}^{(k)} = diag (ã_{ii}^{(k)}) .  (95)

Matrix Ẽ^{(N+1)} obviously consists of elements which have undergone the maximal number of changes. If s(i, j) denotes the number of changes of the element on position (i, j), then

s(i, j) ≤ 2n − 4 , i, j ∈ {1, …, n}, i ≠ j .

The quantity s(i, j) depends upon (i, j) and the pivot strategy. Elements of the matrix Ẽ^{(N+1)} can therefore be denoted as z̃_{s(i,j)}.

Having in mind relation (94), we can now write

| Ẽ^{(N+1)} | ≤ \frac{1}{(d_{N})^{2n−4}} (| P̃^{(2)} || ω̃_{2} | + | P̃^{(3)} || ω̃_{3} | + … + | P̃^{(N)} || ω̃_{N} |) .

(96)

Here the notation | C |= (| c_{ij} |) for C = (c_{ij}) is used. Matrix P̃^{(k)} consists precisely of those elements of l(k)–th and m(k)–th row and column of the matrix Ẽ^{(k)} which already were pivot elements^{3}, i. e. of elements which contribute to the final estimate. All other elements of the matrix P̃^{(k)} are zeros.

Assertion (i) of Lemma 15 gives us

∥ | P̃^{(k)} | ∥=∥ P̃^{(k)} ∥≤ S(Ã^{(k)}) ≤ \sqrt{1.566} · S(Ã^{(1)}) , k = 2, …, N .

(97)

From the relations (96) and (97), Lemma 15 and the Cauchy–Schwarz inequality we obtain

S(Ã^{(N+1)}) =∥ Ẽ^{(N+1)} ∥ = ∥ | Ẽ^{(N+1)} | ∥≤ \frac{1}{(d_{N})^{2n−4}} \sqrt{1.566} · S(Ã^{(1)}) ∑_{k=2}^{N} | ω̃_{k} |

≤ \frac{1.252}{(d_{N})^{2n−4}} · S(Ã^{(1)}) · [(N − 1) ∑_{k=2}^{N} ω̃_{k}^{2}]^{\frac{1}{2}} ≤ \frac{1.252}{(d_{N})^{2n−4}} · S(Ã^{(1)}) · [N ∑_{k=1}^{N} ω̃_{k}^{2}]^{\frac{1}{2}} .(98)

Since N ≥ 3, from Lemma 10 follows that

\frac{1}{(d_{N})^{2n−4}} = \frac{1}{(1 − 0.077\frac{2}{3N})^{n−2}} < 1 + \frac{12}{7} · 0.077 · \frac{2}{3N} (n − 2)

≤ 1 + \frac{12}{7} · 0.077 · \frac{4}{3} · \frac{1}{n} ≤ 1.059 .

Finally, inserting this inequality and assertion (ii) of Lemma 15 into relation (98), we obtain

S(Ã^{(N+1)}) ≤ 0.4 · S(Ã^{(1)})\sqrt{N} · \frac{ε̃_{1}}{δ} .

Applying a similar analysis to matrices B̃^{(k)} yields

S(B̃^{(N+1)}) ≤ 0.4 · S(B̃^{(1)})\sqrt{N} · \frac{ε̃_{1}}{δ} .

From the last two inequalities and the definitions of ε̃_{N+1} and ε̃_{1} follows

ε̃_{N+1} ≤ 0.4 · \sqrt{N} · \frac{ε̃_{1}^{2}}{δ} ,

3Here (l(k), m(k)) denotes pivot pair in the k–th step so this k should not be confused with the k that was fixed at the beginning of the proof.

and the theorem is proved. Q.E.D.

Note that in the proof of Theorem 16 it is not necessary to assume that the affiliation is preserved, i.e. that the pairs [a_{ii}^{(k)}, b_{ii}^{(k)}] approximate the eigenvalues λ_{i} for i = 1, …, n, k = 1, …, N. However, for large enough k this fact follows from Theorem 17.

æ From Theorem 16 and the assumptions (A1) and (A2) follows that

ε̃_{N+1} < \sqrt{N} · \frac{1}{2N} · ε̃_{1} = \frac{1}{2\sqrt{N}} · ε̃_{1} < 0.3 · ε̃_{1} .  (99)

Applying inductively the relation (99) we obtain

ε̃_{rN+1} ≤ (0.3)^{r} · ε̃_{1} ,  r ≥ 1 .  (100)

Therefore,

lim_{r→∞} ε̃_{rN+1} = 0 .  (101)

From the relation (101) and the assertion (i) of Lemma 15 we conclude that

lim_{k→∞} ε̃_{k} = 0 .  (102)

The relation (102) and Theorem 16 imply the quadratic convergence of the Falk–Langemeyer method according to Definition 7 if the eigenvalues are simple and the pivot strategy is cyclic.

Next we prove that under assumptions of Theorem 16 the sequences of matrices (A^{(k)}, k ≥ 1) and (B^{(k)}, k ≥ 1), generated by the Falk–Langemeyer method, converge towards diagonal matrices.

Theorem 17 Let the assumptions of Theorem 16 hold. Then

lim_{k→∞} A^{(k)} = D_{A} ,  lim_{k→∞} B^{(k)} = D_{B} ,

where D_{A} and D_{B} are diagonal matrices.

Proof. The relation (44) implies that

A^{(k)} = (D_{k})^{−1}Ã^{(k)}(D_{k})^{−1} ,  B^{(k)} = (D_{k})^{−1}B̃^{(k)}(D_{k})^{−1} , (103)

where diagonal matrices D_{k} are defined with the relations (45) and (46). It is therefore sufficient to proove that the sequences (Ã^{(k)}, k ≥ 1) , (B̃^{(k)}, k ≥ 1) and ((D_{k})^{−1}, k ≥ 1) converge towards diagonal matrices. The relation (102) implies that the off-diagonal elements of matrices Ã^{(k)} and B̃^{(k)} tend to zero as k → ∞. Therefore, it remains to proove that for i = 1, …, n the sequences (ã_{ii}^{(k)}, k ≥ 1) and (b̃_{ii}^{(k)}, k ≥ 1) converge. The relation (18) and the assumption (A1) imply that for each k ≥ 1 there exists an ordering of the eigenvalues λ_{i} = [s_{i}, c_{i}], i = 1, …, n, such that

| c_{i}ã_{ii}^{(k)} − s_{i}b̃_{ii}^{(k)} | ≤ \frac{ε̃_{k}^{2}}{2δ} ,  i = 1, …, n .  (104)

Let us consider unit vectors [s_{i}, c_{i}]^{T} and [ã_{ii}^{(k)}, b̃_{ii}^{(k)}]^{T} in R^{2}. The left-hand side of the inequality (104) is | sin φ_{i}^{(k)} | where φ_{i}^{(k)} is the angle between these two vectors. The relations (102) and (104) imply

lim_{k→∞} sin φ_{i}^{(k)} = 0 , i =, 1 …, n .

Hence, for each i the sequence of vectors ([ã_{ii}^{(k)}, b̃_{ii}^{(k)}]^{T} , k ≥ 1) has only finite number of accumulation points in R^{2}. Therefore, it suffices to show that for large enough k the changes in ã_{ii}^{(k)} and b̃_{ii}^{(k)} are arbitrary small. From the relation (102) and Lemma 11 we see that α̃_{k} → 0 and β̃_{k} → 0 as k → ∞. Therefore, the changes in ã_{ii}^{(k)} and b̃_{ii}^{(k)} tend to zero as k → ∞. This prooves that for each i ∈ {1, …, n} limits lim_{k→∞} ã_{ii}^{(k)} and lim_{k→∞} b̃_{ii}^{(k)} exist.

We shall now prove that ((D_{k})^{−1}, k ≥ 1) is a convergent sequence. Looking at the definition of D_{k} (relation (45)) we see that it suffices to prove that for each i ∈ {1, …, n} the sequence (d_{i}^{(k)}, k ≥ 1), converges to a nonzero number. From Proposition 9 we have

d_{i}^{(k)} = d_{i}^{(1)}\overline{d}_{i}^{(2)} ⋯ \overline{d}_{i}^{(k)} ,  i = 1, …, n ,  k ≥ 2 .

From the definiteness of pairs (A^{(1)}, B^{(1)}) and (\overline{A}^{(k)}, \overline{B}^{(k)}) we conclude that d_{i}^{(1)} and \overline{d}_{i}^{(k)} are different from zero for all i and k. Hence it suffices to prove that the infinite product ∏_{k=2}^{∞} \overline{d}_{i}^{(k)} converges^{4}. This product converges if and

only if the product ∏_{k=2}^{∞} (\overline{d}_{i}^{(k)})^{4} converges. Therefore, it suffices to show that the later product is absolutely convergent. From the relation (78) we see that

4Since all factors in the product are nonzero the limit, if exists, is also nonzero.

for i ∈ \{1, …, n\} and k ≥ 2 we can write (\overline{d}_{i}^{(k)})^{4} = 1 + u_{i}^{(k)}, so it suffices to show that the series ∑_{k=2}^{∞} u_{i}^{(k)} are absolutely convergent for all i ∈ {1, …, n}.

The relation (83) of of Lemma 12 implies that

(\overline{d}_{i}^{(k+1)})^{4} ≥ 1 − 0.077 · \frac{ε̃_{k}}{δ} ,  1 = 1, …, n ,  k ≥ 1 .

Looking for upper bound instead of lower bound in the relation (78) and making similar estimates as in the relation (83), we obtain

(\overline{d}_{i}^{(k+1)})^{4} ≤ 1 + 0.077 · \frac{ε̃_{k}}{δ} ,  1 = 1, …, n ,  k ≥ 1 .

Therefore,

| u_{i}^{(k+1)} | = | (\overline{d}_{i}^{(k+1)})^{4} − 1 | ≤ 0.077 · \frac{ε̃_{k}}{δ} ,  1 = 1, …, n ,  k ≥ 1 .

Hence it suffices to show that the series ∑_{k=1}^{∞} ε̃_{k} converges. From the assertion (i) of Lemma 15 we have

ε̃_{rN+i} ≤ 1.3 · ε̃_{rN+1} ,  1 ≤ i ≤ N ,  r ≥ 1 ,

hence it suffices to prove the convergence of the sequence ∑_{r=1}^{∞} ε̃_{rN+1}. From the relation (100) we see that the later series is majorized by the convergent series ∑_{r=1}^{∞}(0.3)^{r} · ε̃_{1}. This proves the absolute convergence of the series ∑_{k=2}^{∞} u_{i}^{(k)} for i ∈ {1, …, n} and therefore the convergence of the sequence ((D_{k})^{−1}, k ≥ 1).

Q.E.D.

Note that the global convergence (i.e. the convergence for all definite pairs (A, B) ) of the Falk–Langemeyer method in the case of cyclic pivot strategies is not yet proved.

We end this section by showing that our asymptotic assumptions also imply ultimate quadratic reduction of ε_{rN+1}. Indeed, for r ≥ 1 the relation (103) implies

ε_{rN+1} ≤ (d_{max}^{(rN+1)})^{2} · ε̃_{rN+1} ,  ε̃_{rN+1} ≤ \frac{1}{(d_{min}^{(rN+1)})^{2}} ε_{rN+1} ,

where

d_{max}^{(rN+1)} = max{d_{1}^{(rN+1)}, …, d_{n}^{(rN+1)}} ,
d_{min}^{(rN+1)} = min{d_{1}^{(rN+1)}, …, d_{n}^{(rN+1)}} .

Theorem 16 implies

ε_{(r+1)N+1} ≤ (d_{max}^{((r+1)N+1)})^{2}ε̃_{(r+1)N+1} ≤ (d_{max}^{((r+1)N+1)})^{2} \frac{\sqrt{N}}{δ} ε̃_{rN+1}^{2}
≤ [\frac{d_{max}^{((r+1)N+1)}}{(d_{min}^{(rN+1)})^{2}} ]^{2} \frac{\sqrt{N}}{δ} ε_{rN+1}^{2} ≤ c · \frac{\sqrt{N}}{δ} ε_{rN+1}^{2} ,  r ≥ 1 ,

where c is an upper bound of the convergent sequence ([d_{max}^{((r+1)N+1)}/(d_{min}^{(rN+1)})^{2}]^{2}, r ≥ 1). In a similar way one can prove that quadratic reduction of ε_{rN+1} ultimately implies quadratic reduction of ε̃_{rN+1}.

The techniques described in this section can be used for studying asymptotic convergence properties of various different Jacobi–type algorithms.

5 Concluding remarks

In Algorithm 4 only (l, m)−restrictions of the pair (A^{(k)}, B^{(k)}) are used in each step. Therefore, parallel strategies are in fact cyclic (see [10]) and Theorem 16 and Theorem 17 hold for them as well.

In [13] it is proved that if the assumptions of Theorem 16 hold and the pivot strategy is serial, then

ε̃_{N+1} ≤ \frac{ε̃_{1}^{2}}{δ} .

Modified method. If the problem (1) has multiple eigenvalues, the method can fail to be quadratically convergent. This failure occurs because when pairs [a_{ll}^{(k)}, b_{ll}^{(k)}] and [a_{mm}^{(k)}, b_{mm}^{(k)}] (here (l, m) is the pivot pair in the k–th step) approximate the same eigenvalue, then parameters α̃_{k} and β̃_{k} can be of order O(1) and, therefore, some previously annihilated elements can become of order O(ε̃_{k}) again. This situation is described in detail in [7] and [13]. Simple omitting of these critical steps does not yield to the quadratic convergence, even though the measure τ̃_{k} = τ(Ã^{(k)}, B̃^{(k)}), k ≥ 1 , from Corollary 2 tends to zero. The relation (16) does not imply that the off-diagonal elements of diagonal blocks tend to zero together with τ̃_{k}, but merely that the diagonal blocks become more and more proportional. Therefore, ε̃_{k} does not have to tend to zero at all and the convergence of τ̃_{k} can considerably slow down. If we modify the method so that in such cases we use triangular transformation matrices similar to the matrix from step (5a) of Algorithm 4, the quadratic convergence persists.

Modification of the Falk–Langemeyer method and the proof of quadratic convergence of the modified method will be topics of our subsequent paper.

Numerical results. Our test program is written in FORTRAN in double precision. Test pairs were generated in the manner that A = G^{T} D_{A}G and B = G^{T} D_{B}G, where diagonal matrices D_{A} and D_{B} are being read and G is random. For elements of matrix G only numbers which are sums of the powers of 2 were used, so the test pairs were stored as accurately as possible.

The iterative process is terminated when, after some cycle r, inequality

ε_{rN+1} < eps · \sqrt{∥A∥^{2} + ∥B∥^{2}} · 2N

is fulfilled, where eps is machine precision. After the end of the process, the maximal error of the residual

max_{1≤i≤n} \{\frac{∥ b′_{i}Af′_{i} − a′_{i}Bf′_{i} ∥_{max}}{\sqrt{(a′_{i})^{2} + (b′_{i})^{2}}\sqrt{∥ Af′_{i} ∥^{2} + ∥ Bf′_{i} ∥^{2}}} \} ,

is calculated. Here [a′_{i}, b′_{i}] are the calculated eigenvalues of the pair (A, B) and f′_{i} are the corresponding eigenvectors. Also the maximal absolute values of the off-diagonal elements of matrices (F′)^{T} AF′ and (F′)^{T} BF′ are calculated. Those three quantities were usually of order. Infinite eigenvalues were represented with numbers of order of magnitude O(1/machine precision).

We observed the convergence of both measures ε_{k} and ε̃_{k}. Observations confirmed all theoretical results. For starting pairs that were not almost diagonal, convergence was in the beginning linear and several cycles were needed before quadratic convergence started. The asymptotic assumption (A2) appears to be very adequate because in almost all cases quadratic convergence started after it was fulfilled. Algorithm behaved very regularly in the sense that the condition ℑ^{(k)} ≥ 0, k ≥ 1, (see assertion (i) of Proposition 3) was always fulfilled for definite starting pairs. This condition was fulfilled even in some cases when the starting pair was semidefinite, or slightly indefinite.

Average number of cycles for smaller matrices (n ≤ 15) was around ten and for larger matrices (n ≤ 100) around fifteen. Last cycles were usually empty, i.e. not all N steps were executed. For orientation, the approximate duration of the process is five minutes for n = 40 and one and a half hour for n = 100 on IBM PC/AT with a coprocessor, and about thirty times shorter on IBM 4371.

In the presence of very close eigenvalues several additional cycles were usually needed because the quadratic convergence was delayed. The existance of additional cycles does not disagree with theoretical results since the quantity δ from the asymptotic assumption (A2) is in this case very small.

We observed that the results are generally better if increasing or decreasing order of numbers defined with diagonal pairs [a_{ii}^{(k)}, b_{ii}^{(k)}] is preserved by interchanging pivot rows and columns if necessary. However, interchanging must be stopped after the asymptotic assumption (A2) is fulfilled. Otherwise some off–diagonal element which was not yet annihilated can “run away” from annihilation and therefore terminate quadratic convergence.

æ Example. We give an example of the pair of order 10 generated in the previously described manner. Elements of the matrices D_{A} and D_{B} are

−2, 1, 10, 0, −0.001, 10, 1, 5, 5, 4

and

−1, 0.1, −1, −100, −100, 0, −1, 0.1, 1, 1,

respectively, so the exact eigenvalues of the problem are

2, 10, −10, 0, 0.00001, ∞, −1, 50, 5, 4 .

Elements of the matrix G are uniformly distributed integers from the interval [−10, 10]. Note that both matrices A and B are indefinite, while the pair (A, B) itself is definite (for example A − 3B > 0). In order to increase the stability of the computation, the process started from the normalized pair (Ã, B̃).

Only upper triangles of the matrices A and B are displayed. Each row begins with the diagonal element. Asymptotic convergence is described as follows: in column CYC is the number of cycle; in column ROT is the number of rotations performed in the cycle; columns SUMA, SUMB, SUM and SUMT display values of S(A^{(k)}), S(B^{(k)}), ε_{k} and ε̃_{k} after the cycle, respectively.

ORDER OF MATRICES N =  10
COLUMN CYCLIC PIVOT STRATEGY
STOPPING CRITERION: SUM(K) <  .49D-13
MATRIX A
ROW
1    .21350D+04    .41900D+03    .11600D+03   -.11430D+04   -.10490D+04
     .44002D+03   -.13750D+04    .20027D+02    .51903D+03   -.60802D+03
2    .14310D+04   -.32700D+03   -.34200D+03   -.26100D+03   -.10390D+04
    -.29600D+03   -.43200D+03    .50000D+03   -.13100D+03
3    .18320D+04    .28000D+03    .93300D+03    .64100D+03   -.91000D+03
    -.38500D+03   -.74500D+03   -.58200D+03
4    .11860D+04    .85799D+03    .34099D+03    .56100D+03   -.40001D+03
    -.49101D+03    .57401D+03
5    .99295D+03    .68695D+03    .22402D+03   -.53006D+03   -.69606D+03
    -.90965D+02
6    .13360D+04   -.57298D+03   -.43106D+03   -.97606D+03    .43035D+02
7    .13470D+04    .29703D+03    .10027D+02    .67399D+03
8    .87292D+03    .32992D+03   -.49550D+01
9    .88792D+03    .40105D+03
10    .82798D+03

MATRIX B
ROW
1   -.72420D+04   -.81550D+04    .40130D+04   -.20630D+04   -.39800D+03
     .31460D+04   -.73170D+04    .83650D+04    .27080D+04   -.78930D+04
2   -.99425D+04    .50853D+04   -.28814D+04   -.28007D+04    .93320D+03
    -.79655D+04    .68774D+04   -.11530D+03   -.79841D+04
3   -.26020D+04    .14878D+04    .14297D+04   -.63790D+03    .39712D+04
    -.34409D+04   -.16500D+02    .39654D+04
4   -.10848D+04   -.15558D+04   -.55260D+03   -.20896D+04    .11306D+04
    -.83980D+03   -.18888D+04
5   -.59091D+04   -.46836D+04   -.23130D+03   -.41948D+04   -.61837D+04
     .10755D+04
6   -.50607D+04    .28013D+04   -.69463D+04   -.62826D+04    .43479D+04
7   -.72920D+04    .82531D+04    .26803D+04   -.78484D+04
8   -.12938D+05   -.81476D+04    .10050D+05
9   -.82005D+04    .46091D+04
10   -.88056D+04

A S Y M P T O T I C   C O N V E R G E N C E
 CYC   ROT    SUMA     SUMB     SUM      SUMT
  1     45   .60D+00  .35D+01  .35D+01  .18D+01
  2     45   .72D+00  .32D+00  .79D+00  .12D+01
  3     45   .57D+00  .28D+00  .64D+00  .72D+00
  4     45   .31D+00  .21D+00  .38D+00  .24D+00
  5     45   .18D-01  .23D-01  .30D-01  .25D-01
  6     45   .39D-02  .13D-01  .13D-01  .93D-02
  7     45   .56D-04  .15D-03  .16D-03  .92D-04
  8     44   .18D-09  .61D-09  .63D-09  .19D-09
  9     29   .26D-20  .19D-20  .32D-20  .93D-20
TOTAL NO. OF ROTATIONS   388          TIME(sec)    5.68

C A L C U L A T E D   E I G E N V A L U E S
I       A(I,I)            B(I)                D(I)
1    .36774301D+01    .91935754D+00    .40000000000001D+01
2    .29605919D+01    .59211837D-01    .50000000000001D+02
3    .24951699D+00    .24951699D-01    .10000000000015D+02
4    .62005903D+00    .12401181D+00    .50000000000011D+01
5   -.16044487D+00   -.80222433D-01    .19999999999999D+01
6   -.94464581D-04   -.94464581D+01    .10000000000004D-04
7   -.43510239D-16   -.23229057D+01    .18730953468577D-16
8    .62365226D-01   -.62365226D-01   -.99999999999990D+00
9    .42815492D+01   -.25722206D-13   -.16645341957511D+15
10    .32339485D+01   -.32339485D+00   -.99999999999999D+01
 MAXIMAL(relative) ERROR =  .15D-13   FOR  I =   3
 MAXIMAL OFF-DIAGONAL ELEMENTS:
 Ft A F =  .35D-14        Ft B F =  .30D-13

References

[1] Charlier, J.-P., Van Dooren, P., A Jacobi–like algorithm for computing the generalized Schur form of a regular pencil, Philips Research Laboratory, technical report,1988.

[2] Falk, S., Langemeyer, P., Das Jacobische Rotations–Verfahren für reelsymmetrische Matrizen–Paare I, II, Elektronische Datenverarbeitung (1960) 30–43.

[3] Gose, G., Das Jacobi Verfahren für Ax = λBx , ZAMM 59(1979) 93–101.

[4] Hari, V., On Cyclic Jacobi Methods for the Positive Definite Generalized Eigenvalue Problem, Disertation, Fernuniversität Hagen, 1984.

[5] Hari, V., On the Convergence of cyclic Jacobi–Like Processes, Lin. Alg. and Its Appl. 81(1986) 105–127.

[6] Hari, V., On the Quadratic Convergence of Jacobi Algorithms, Radovi matematički 2(1986) 127–146.

[7] Hari, V., On Pairs of Almost Diagonal Matrices, proposed for publication in Lin. Alg. and Its Appl.

[8] Henrici, P., Zimmermann, K., An Estimate for the Norms of Certain Cyclic Jacobi Operators, Lin. Alg. and Its Appl. 1(1968) 489–501.

[9] Kempen, H.P.M. van, On the Quadratic Convergence of the Serial Cyclic Jacobi Method, Numer. Math. 9(1966) 19–22.

[10] Luk, F., Park, H., On the equivalence and convergence of parallel Jacobi algorithms, Proceedings SPIE conference on Advanced Algorithms and Architecture II, 1987.

[11] Parlett, B.N., Symmetric Eigenvalue Problem, Prentice Hall Inc., Englewood Cliffs, N.J., 1980.

[12] Sameh, A. H., On Jacobi and Jacobi–like Algorithms for a Parallel Computer, Mathematics of Computation, Vol. 25 No. 9(1971) 579–590.

[13] Slapničar, I., Quadratic Convergence of the Falk–Langemeyer Method, M. S. Thesis, University of Zagreb, 1988. (in Croatian language)

[14] Stewart, G.W., Perturbation Bounds for the Definite Generalized Eigenvalue Problem, Lin. Alg. and Its Appl. 23(1960) 69–85.

[15] Veselić, K., An eigenreduction algorithm for definite matrix pairs and its applications to overdamped linear systems, to appear in SIAM J. Sci. Stat. Comp.

[16] Wilkinson, J.H., Note on the Quadratic Convergence of the Cyclic Jacobi Processes, Numer. Math. 4(1962) 296–300.

[17] Wilkinson, J.H., Almost Diagonal Matrices with Multiple or Close Eigenvalues, Lin. Alg. and Its Appl. 1(1968) 1–12.

[18] Wilkinson, J.H., The Algebraic Eigenvalue Problem, Oxford University Press, Oxford,1965.

[19] Zimmermann, K., On the Convergence of a Jacobi Process for Ordinary and Generalized Eigenvalue Problems, Disertation No 4305(1965), Eidgenossische Technische Hochschule Zürich.

æ