On the Quadratic Convergence of the Paardekooper Method II

Vjeran Hari

Department of Mathematics, University of Zagreb,

Bijenička cesta 30, 41001 Zagreb, Croatia and Noah H. Rhee Department of Mathematics, University of Missouri-Kansas City, Kansas City, MO 64110-2499

Abstract

In this paper we consider the asymptotic convergence of the Paardekooper method for reducing a skew-symmetric matrix to its Murnaghan form. Using a new technique we derive a sharp constant bound for the quadratic convergence of the method under the rowcyclic strategy. We provide numerical examples which show that the proposed asymptotic conditions are not only sufficient but in the general case also necessary.

Subject Classifications: AMS (MOS): 65F15; CR: 5.14

0 Introduction

With prevailing trend of parallel processing diagonalization methods for solving eigenvalue problems came into renewed interest of numerical analysts. This interest is due to an inherent parallelism of such methods and also to their ability to deliver the computed eigenvalues and eigenvectors with high relative accuracy [1]. Here we study a method for reducing a skew-symmetric matrix to its Murnaghan form [9]. The special feature of the method is that it operates on 2 by 2 blocks and thus is more efficient than the standard symmetric Jacobi method. It is especially efficient on the machines for which simple 4 by 4 matrix transformation instructions are available (e.g. IIT coprocessor support software). The method is due to M.H.C. Paardekooper [10]. We note that skew-symmetric eigenvalue problems appear in different areas of applied mathematics, e.g. structural mechanics [2, 12].

Incentives for writing this paper came with discovery of a new, more powerful technique for proving asymptotic convergence bounds of Jacobi-type methods [6, 11]. We wanted to show that the new technique can successfully be applied to a method which essentially works with blocks, what usually complicates the analysis and worsens the estimates.

The existing result on Paardekooper method [7] considers the quadratic convergence under a general cyclic strategy in the case of distinct eigenvalues. It is based on the well known Wilkinson’s paper [13] for symmetric matrices. An attempt [5] to make a similar analysis for the row-cyclic method brought a bound depending on n, the matrix order, and the proof is made for simple eigenvalues only.

In this paper we derive a sharp constant quadratic convergence bound of the method under the row-cyclic strategy. The analysis is much more subtle than the one concerning the asymptotic convergence of the symmetric Jacobi method. Our asymptotic assumptions are essentially weaker than the one used in [7, 5] since they allow the case of double eigenvalues or pairs of simple eigenvalues which form clusters. By numerical examples we show that in the case of eigenvalues with multiplicities larger than two or clusters with more than two eigenvalues we can have a failure of the quadratic convergence. The numerical examples also show how very accurate the obtained bound is.

The paper is organized as follows. In Section 1 we define the cyclic Paardekooper method. In Section 2 we derive the basic bounds for the angles under a suitable asymptotic assumption. In Section 3 we derive a sharp quadratic convergence bound for the row-cyclic method. At the end of section 3 we present numerical tests.

1 The Description of the Method

In this section we briefly describe the method and introduce notations. Let A = (a_{ij}) be a skew-symmetric matrix of order n. Using notation from [10] we assume the partition

A = [ \begin{matrix}A_{11} & ⋯ & A_{1k} \\ ⋮ & ... & ⋮ \\ A_{k1} & ⋯ & A_{kk} \end{matrix} ] ,  (1.1)

where each A_{ij} is 2 × 2 except for the edge submatrices A_{ik} (A_{ki}), i = 1, 2, …, k − 1, which are 2 × 1 (1 × 2) when n is odd. Of course, in such a case A_{kk} = [0]. Thus k = [(n + 1)/2] where [a] denotes the largest integer ≤ a. Instead of (1.1) we will also write A = (A_{ij}). With the partition (1.1) is associated

S(A) = [∑_{i<j} ∥A_{ij}∥^{2}]^{1/2} ,  (1.2)

which is called (see [10]) departure from the Murnaghan form. In (1.2) as well as throughout the whole paper ∥ · ∥ stands for the Euclidean norm, that is, ∥A∥^{2} = trace(A^{T} A). Paardekooper method is a Jacobi-type method which reduces A to its Murnaghan form. Hence one step of the method uses an orthogonal similarity transformation

A′ = ℑ_{ij}^{T} A ℑ_{ij}, i < j,

to make A′_{ij} zero. The matrix ℑ_{ij} is called Jacobi annihilator of A_{ij}, and A_{ij} is called pivot submatrix.

If A_{ij} is 2 × 2 then

ℑ_{ij} = R_{ij}^{1}R_{ij}^{2}R_{ij}^{3}R_{ij}^{4}, i < j,

where R_{ij}^{t}, 1 ≤ t ≤ 4 are rotation matrices. If R(p, q; ϕ) stands for the rotation matrix in the (p, q) plane by the angle ϕ, with the four essential elements r_{pp} = r_{qq} = cos ϕ, r_{pq} = −r_{qp} = sin ϕ, then (cf. [10])

R_{ij}^{1} = R(2i − 1, 2j − 1; ϕ_{1}),
R_{ij}^{2} = R(2i, 2j; ϕ_{2}),
R_{ij}^{3} = R(2i − 1, 2j; ϕ_{3}),
R_{ij}^{4} = R(2i, 2j − 1; ϕ_{4}).

The rotation R_{ij}^{1} [R_{ij}^{3}] modifies the elements at positions (2i−1, 2j), (2i,2j − 1) [(2i−1, 2j−1), (2i, 2j)] in such a way that R_{ij}^{2} [R_{ij}^{4}] can annihilate them. The skew-simmetry property ensures that zeros obtained at positions (2i − 1, 2j), (2i, 2j − 1) are not violated by the transformations R_{ij}^{3} and R_{ij}^{4}. We also write

ℑ_{ij} = ℑ_{ij}^{1} ℑ_{ij}^{2}, ℑ_{ij}^{1} = R_{ij}^{1}R_{ij}^{2}, ℑ_{ij}^{2} = R_{ij}^{3}R_{ij}^{4}.

If n = 2k − 1, then A_{ik} is 2 × 1. In this case we set

ℑ_{ik} = R_{ik}^{1}R_{ik}^{2}, 1 ≤ i ≤ k − 1,

with

R_{ik}^{1} = R(2i, n; ϕ_{1}),
R_{ik}^{2} = R(2i − 1, n; ϕ_{2}).

Here R_{ik}^{1} (R_{ik}^{2}) makes the element a_{2i−1,n} (a_{2i,n}) zero.

Since A_{ij} is zeroed and ℑ_{ij} is orthogonal with special stucture we have (see [10])

S^{2}(A′) = S^{2}(A) − ∥A_{ij}∥^{2}.  (1.3)

More precisely, if A_{ij} ∈ R^{2×2}, then

S^{2} ((ℑ_{ij}^{1})^{T} Aℑ_{ij}^{1}) = S^{2}(A) − (v^{2} + w^{2}),  (1.4)
S^{2}(A′) = S^{2} ((ℑ_{ij}^{1})^{T} Aℑ_{ij}^{1}) − (x^{2} + y^{2}).  (1.5)

If A_{ik} ∈ R^{2×1}, then

S^{2} ((R_{ik}^{1})^{T} AR_{ik}^{1}) = S^{2}(A) − v^{2},  (1.6)
S^{2}(A′) = S^{2} ((R_{ik}^{1})^{T} AR_{ik}^{1}) − w^{2}.  (1.7)

The angle formulas can be found in [10, 7].

Besides angle determination we have to specify a pivot strategy. A pivot strategy is a function from N_{0} = {0, 1, …} to P_{k} = {(i, j); 1 ≤ i < j ≤ k}. In this paper We consider mainly the cyclic strategies. They can be defined by the orderings of the set P_{k} = {(i, j); 1 ≤ i < j ≤ k}. If

O = (i_{0}, j_{0}), (i_{1}, j_{1}), …, (i_{K−1}, j_{K−1}),

where K = k(k − 1)/2, which is the number of elements in P_{k}, the appropriate cyclic strategy is defined by the index functions i(r), j(r) as follows: (i(r), j(r)) = (i_{s}, j_{s}) provided that r ≡ s (modK), 0 ≤ s ≤ K −1, r ≥ 0. The most common cyclic strategies are the row- and column-cyclic ones defined by the orderings

O_{R} = (1, 2), (1, 3), …, (1, k); (2, 3), …, (2, k); …; (k − 1, k)

and

O_{C} = (1, 2); (1, 3), (2, 3); …; (1, k), (2, k), …, (k − 1, k),

respectively. The first K applications of Jacobi annihilator ℑ_{ij} according to the given ordering O make the first cycle or sweep of the cyclic Paardekooper method. In this paper we prove a sharp quadratic convergence bound for the row-cyclic Paardekooper method. Now we introduce notations.

Notations

For a given Paardekooper method we use the following notations:

A^{(0)} = A = A^{[0]}: initial skew-symmetric matrix of order n

A^{[r]} = (A_{ij}^{[r]}) = (a_{st}^{[r]}): matrix obtained \underline{after} the r’th step

A^{(r)} = (A_{ij}^{(r)}) = (a_{st}^{(r)}): matrix obtained \underline{after} the r’th rotational transformation. Thus A^{[r]} = A^{(4r)} provided that all annihilated submatrices were 2 × 2.

If only one cycle is considered, then we use the following notations:

A^{[i,j]} = (A_{lm}^{[i,j]}) = (a_{st}^{[i,j]}) : matrix obtained just \underline{before} the transformation with ℑ_{ij}

A^{(p,q)} = (A_{lm}^{(p,q)}) = (a_{st}^{(p,q)}): matrix obtained just \underline{before} applying the rotation matrix in the (p, q) plane

A^{[i,j]+t} = (A_{lm}^{[i,j]+t}) = (a_{st}^{[i,j]+t}): matrix obtained from A^{[i,j]}, \underline{after} performing t additional rotational transformations

A^{(i,j)+t} = (A_{lm}^{(i,j)+t}) = (a_{st}^{(i,j)+t}): matrix obtained from A^{(i,j)}, \underline{after} performing t additional rotational transformations

Block(i, j): = \{ \begin{matrix}\{(2i − 1, n), (2i, n)\} if j = k, n = 2k − 1 \\ \{(2i − 1, 2j − 1), (2i, 2j − 1), (2i − 1, 2j), (2i, 2j)\} otherwise \end{matrix}

ϕ_{pq}: the angle used in the rotation matrix in the (p, q) plane.

Φ_{ij}:  Φ_{ij} = ∑_{(p,q)∈Block(i,j)} sin^{2} ϕ_{pq} .

2 Asymptotic Assumption and Angle Estimates

We denote the eigenvalues of A by ±ıν_{j}, where ν_{j} ≥ 0 for 1 ≤ j ≤ k. Let

3δ = min_{i≠j} |ν_{i} − ν_{j}| > 0.  (2.1)

That is, we are assuming distinct eigenvalues (0 is the only possible double eigenvalue of A). Note that 3δ can be significantly larger than min{|λ_{i} − λ_{j}|; λ_{i} ≠ λ_{j}, λ_{i}, λ_{j} ∈ σ(A)}, where σ(A) is the spectrum of A (see Figure 1 below). Later we shall modify δ to deal with pairs of identical or close eigenvalues.

Asymptotic Assumption

[figure: ν_{4}; ν_{3}; ν_{2}; ^{3δ}; ν_{1}; 0]

Figure 1: Eigenvalues of ıA

We make the following asymptotic assumption.

S(A) < \frac{1}{2}δ.  (2.2)

Note that this assumption is equivalent to the assumption made in [7].

Applying Hoffman–Wielandt theorem (see [8, 7]) to A^{[r]} and diag(A_{11}^{[r]}, …, A_{kk}^{[r]}) from the partition (1.1) and using the relations (1.3) and (2.2) we obtain

∑_{i=1}^{k′} ||a_{2i−1,2i}^{[r]}| − ν_{i}^{r}|^{2} ≤ S^{2}(A^{[r]}) ≤ S^{2}(A) < \frac{δ^{2}}{4} , r ≥ 0,  (2.3)

where k′ = [n/2], and ν_{1}^{r}, …, ν_{k′}^{r} is an appropriate ordering of (positive) ν_{j} if n is (odd) even. Hence

||a_{2i−1,2i}^{[r]}| − ν_{i}^{r}| < \frac{δ}{2}, i = 1, …, k′, r ≥ 0.  (2.4)

We say that 2 × 2 block A_{ii}^{[r]} is affiliated with the eigenpair ±ıν_{j} (or simply with ν_{j}) if ||a_{2i−1,2i}^{[r]}|−ν_{j}| < δ/2. So the relation (2.4) says that A_{ii}^{[r]} is affliated with ν_{i}^{r}, i = 1, …, k′, r ≥ 0. We note that A_{ii}^{[r]} is affliated with none other eigenvalues than ν_{i}^{r}. This can be seen as follows. If ν_{j}^{r} ≠ ν_{i}^{r}, by using the relations (2.1) and (2.4) we obtain

||a_{2i−1,2i}^{[r]}| − ν_{j}^{r}| = ||a_{2i−1,2i}^{[r]}| − ν_{i}^{r} + ν_{i}^{r} − ν_{j}^{r}|
≥ |ν_{i}^{r} − ν_{j}^{r}| − ||a_{2i−1,2i}^{[r]}| − ν_{i}^{r}|
≥ 3δ − \frac{δ}{2} ≥ \frac{5}{2}δ,

which shows that A_{ii}^{[r]} is not affiliated with ν_{j}^{r}.

It turns out that the affliation of the diagonal blocks is preserved in each step of the Paardekooper method under any pivot strategy with the assumption (2.2). So the relation (2.3) can be simplified

∑_{i=1}^{k′} ||a_{2i−1,2i}^{[r]}| − ν_{i}|^{2} ≤ S^{2}(A^{[r]}) ≤ S^{2}(A) < \frac{δ^{2}}{4} , r ≥ 0  (2.5)

The following Lemma justifies it. In the Lemma we use the notations

A_{ii} = [ \begin{matrix}0 & α \\ −α & 0 \end{matrix} ] , A_{ij} = [ \begin{matrix}x & v \\ w & y \end{matrix} ] , A_{jj} = [ \begin{matrix}0 & β \\ −β & 0 \end{matrix} ] .  (2.6)

and

A_{ii}^{(2)} = [ \begin{matrix}0 & α′ \\ −α′ & 0 \end{matrix} ] , A_{jj}^{(2)} = [ \begin{matrix}0 & β′ \\ −β′ & 0 \end{matrix} ] .  (2.7)

Lemma 2.1 If the assumption (2.2) holds, the affiliation of the diagonal blocks is preserved in each step under any pivot strategy.

Proof. Because of the relation (1.3) it is enough to consider only the first step. Let A_{ii}, A_{jj} and A_{ij}, i < j be the essential blocks. We assume that A_{ii} (A_{jj}) is affiliated with ν_{i} (ν_{j}) (see the relation (2.5) with r = 0), that is, we are assuming

||α| − ν_{i}| < \frac{δ}{2}, ||β| − ν_{j}| < \frac{δ}{2}.

First we consider the case when A_{jj} is of order two. Since the relations (1.4) and (1.5) imply S(A^{[1]}) ≤ S(A^{(2)}) ≤ S(A) we can consider the transformations with ℑ_{ij}^{1} and ℑ_{ij}^{2} separately. Both transformations are treated in the same way. So it is enough to consider only the transformation with ℑ_{ij}^{1}. Using the condition

(c_{1}v − s_{1}β)(c_{1}β + s_{1}v) = (c_{1}w − s_{1}α)(c_{1}α + s_{1}w), c_{1} = cos ϕ_{1}, s_{1} = sin ϕ_{1}

which follows from the requirement a_{2i−1,2j}^{(2)} = 0 = a_{2i,2j−1}^{(2)}, we obtain

α′ = (c_{1}α + s_{1}w)c_{2} − (c_{1}v − s_{1}β)s_{2}, c_{2} = cos ϕ_{2}, s_{2} = sin ϕ_{2}
= (c_{1}α + s_{1}w)c_{2}(1 + t_{2}^{2}),  t_{2} = tan ϕ_{2}
= \frac{c_{1}}{c_{2}} α + \frac{s_{1}}{c_{2}} v,

β′ = (c_{1}β + s_{1}v)c_{2} − (c_{1}w − s_{1}α)s_{2},
= (c_{1}β + s_{1}v)c_{2}(1 + t_{2}^{2}),
= \frac{c_{1}}{c_{2}} β + \frac{s_{1}}{c_{2}} v.

Since all other diagonal blocks are invariant under the transformation ℑ_{ij}^{1}, it follows that affiliation is preserved provided that

(α − β)(α′ − β′) > 0.

Since

α′ − β′ = \frac{c_{1}}{c_{2}} (α − β) − \frac{s_{1}}{c_{2}} (w − v)

we have

(α − β)(α′ − β′) = \frac{c_{1}}{c_{2}} (α − β)^{2} + \frac{s_{1}}{c_{2}} (α − β)(w − v)
≥ \frac{c_{1}}{c_{2}} |α − β|(|α − β| − |w − v||t_{1}|), t_{1} = tan ϕ_{1}
> \frac{c_{1}}{c_{2}} 2.29 δ ( 2.29δ − \sqrt{2}∥A_{ij}∥ · 1)
> \frac{c_{1}}{c_{2}} 2.29 δ ( 2.29 − \frac{\sqrt{2}}{2})δ > 0

Here we used the following inequality.

|α − β| ≥ ||α| − |β||
≥ |ν_{i} − ν_{j}| − ||α| − ν_{i}| − ||β| − ν_{j}|
> 3 δ − \sqrt{2}\sqrt{(|α| − ν_{i})^{2} + (|β| − ν_{j})^{2}}
≥ (3 − \frac{\sqrt{2}}{2}) δ > 2.29289δ  (2.8)

If A_{ik} is 2 × 1, then A_{kk} = [0]. So, A_{ii}^{[1]} is necessarily affiliated with ν_{i}. Otherwise, we would get a contradiction to the relation (2.3).

Angle Estimates

Here we derive bounds for the angles ϕ_{t}, t = 1, 2, 3, 4 of a Jacobi annihilator ℑ_{ij}. In the case when n is odd and j = k we estimate the angles ϕ_{1} and ϕ_{2}. In the following Lemma again we use the notations defined in the

relations (2.6) and (2.7) together with A_{ik} = [ \begin{matrix}v \\ w \end{matrix} ] when n is odd.

Lemma 2.2 Let the condition (2.2) hold.

(i) If A_{ij} is 2 × 2 then

| sin ϕ_{t} |≤ \{ \begin{matrix}c · \frac{\sqrt{v^{2}+w^{2}}}{δ} & , & t = 1 \\ b · \frac{\sqrt{v^{2}+w^{2}}}{δ} & , & t = 2 \\ c · ^{\frac{\sqrt{x^{2}+y^{2}}}{δ}} & , & t = 3 \\ b · ^{\frac{\sqrt{x^{2}+y^{2}}}{δ}} & , & t = 4 \end{matrix}  (2.9)

and

Φ_{ij} ≤ a^{2} · \frac{∥A_{ij}∥^{2}}{δ^{2}}

where

a = 1.0062, b = 0.89693, c = 0.45602.

(ii) If n = 2k − 1 and j = k (in this case A_{ik} is 2 × 1) then

| sin ϕ_{t} |≤ \{ \begin{matrix}c′ · \frac{|v|}{δ} & , & t = 1 \\ b′ · \frac{|w|}{δ} & , & t = 2 \end{matrix}  (2.10)

and

Φ_{ik} ≤ (a′)^{2} · \frac{∥A_{ik}∥^{2}}{δ^{2}}  (2.11)

where

a′ = 0.425, b′ = 0.425, c′ = 0.4

Proof. (i) From the angle formula [7] we obtain

| sin ϕ_{1}| ≤ | tan ϕ_{1}| ≤ \frac{1}{2}| tan 2ϕ_{1}| ≤ \frac{|α||w| + |β||v|}{|α^{2} − β^{2}| − |v^{2} − w^{2}|}
≤ \frac{\sqrt{v^{2} + w^{2}}}{\frac{|α|+|β|}{\sqrt{α^{2}+β^{2}}} ||α| − |β|| − \frac{max\{v^{2},w^{2}\}}{max\{|α|,|β|\}}}

≤ \frac{\sqrt{v^{2} + w^{2}}}{||α| − |β|| − \frac{∥A_{ij}∥^{2}}{max\{|α|,|β|\}}} .

Since A_{ii} (A_{jj}) is affiliated with ν_{i} (ν_{j}), by suing the relation (2.1) we obtain

max\{|α|, |β|\} ≥ \{ \begin{matrix}ν_{i} − |ν_{i} − |α|| & , & ν_{i} > ν_{j} \\ ν_{j} − |ν_{j} − |β|| & , & ν_{i} < ν_{j} \end{matrix} \}
≥ 3 δ − \frac{1}{2}δ ≥ \frac{5}{2}δ,  (2.12)

From the relations (2.12), (2.8) and (2.2) we obtain

||α| − |β|| − \frac{∥A_{ij}∥^{2}}{max\{|α|, |β|\}} >

2.29289 δ − \frac{S^{2}(A)}{2.5δ} > (2.29289 − 0.1) δ > 2.19289δ

and since 1/2.19289 < 0.45602, the constant c from the relation (2.9) is proved for t = 1.

In order to obtain a bound for sin ϕ_{2} we consider two cases: |α| ≥ |β| and |α| < |β|. In the first (second) case we use the first (second) formula for tan ϕ_{2} in [7]. We shall consider only the case |α| ≥ |β|. The second case is treated similarly.

We have

| sin ϕ_{2}| ≤ | tan ϕ_{2}| ≤ \frac{|v| + |β|| tan ϕ_{1}|}{|α| − |w|| tan ϕ_{1}|}

≤ \frac{\frac{|v|}{|α|} + \frac{|β|}{|α|}| tan ϕ_{1}|}{1 − \frac{|w|}{|α|}| tan ϕ_{1}|} ≤ \frac{\frac{|v|}{2.5δ} + 1 · 0.45602 · \frac{\sqrt{v^{2}+w^{2}}}{δ}}{1 − \frac{|w|}{2.5δ} · 0.45602\frac{\sqrt{v^{2}+w^{2}}}{δ}}

≤ \frac{(\frac{1}{2.5} + 0.45602) \frac{\sqrt{v^{2}+w^{2}}}{δ}}{1 − \frac{0.45602}{2.5} \frac{1}{4}}
≤ 0.89693\frac{\sqrt{v^{2} + w^{2}}}{δ} .

After the first two rotational transformations we see from (1.4) that the relations (2.8) and (2.12) still hold. Note the formulas for tan 2ϕ_{3} and tan ϕ_{4} (see [7]) have the same form as those for tan 2ϕ_{1} and tan ϕ_{2}, respectively. Besides the proof of Lemma 2.1 shows the affiliation of the diagonal blocks have not changed after the first two rotational transformations. So we can use exactly the same proof as above to obtain the remaining two inequalities in the relation (2.9). All we need to change is to replace α, β, v and w by α′, β′, x and y, respectively.

Squaring and summing the inequalities in (2.9) and using the notation in the relation (2.6) we obtain

Φ_{ij} ≤ (0.45602^{2} + 0.89693^{2})\frac{∥A_{ij}∥^{2}}{δ^{2}} ≤ 1.0062^{2} \frac{∥A_{ij}∥^{2}}{δ^{2}} .

(ii) The proof of this part is given in [7].

Note that the relation (2.9) is a refinement of the results in [7]. We also note that in any case

Φ_{ij} ≤ 1.0062^{2} \frac{∥A_{ij}∥^{2}}{δ^{2}} ,

since a′ = 0.425 < 1.0062 = a.

Lemma 2.3 Let A be a skew-symmetric matrix satisfying (2.2). Let the sequence A^{[0]} = A, A^{[1]}, …, A^{[K]} be obtained by applying Paardekooper method to A under any cyclic strategy. Let Φ_{i(r),j(r)} be associated with the pivot submatrix A_{i(r),j(r)}^{[r]}, 0 ≤ r ≤ K − 1. Then

∑_{r=0}^{K−1} Φ_{i(r),j(r)} ≤ [1.0062\frac{S(A)}{δ}]^{2} .  (2.13)

Proof. Since the affiliation of diagonal blocks does not change in each step under any cyclic strategy, A_{i(r),i(r)}^{[r]} and A_{j(r),j(r)}^{[r]} are affiliated with ν_{i(r)} and ν_{j(r)}, respectively. Hence by following the proof of Lemma 2.2 we see that

Φ_{i(r),j(r)} ≤ 1.0062^{2} \frac{∥A_{i(r),j(r)}^{[r]}∥^{2}}{δ^{2}} . 0 ≤ r ≤ K − 1

In fact, all we need to modify in the proof of Lemma 2.2 is to use ν_{i(r)} and ν_{j(r)} in place of ν_{i} and ν_{j}, repectively, and to use the following notations.

A_{i(r),i(r)}^{[r]} = [ \begin{matrix}0 & α \\ −α & 0 \end{matrix} ] , A_{i(r),j(r)}^{[r]} = [ \begin{matrix}x & v \\ w & y \end{matrix} ] , A_{j(r),j(r)}^{[r]} = [ \begin{matrix}0 & β \\ −β & 0 \end{matrix} ] ,

and

A_{i(r),i(r)}^{[r]+2} = [ \begin{matrix}0 & α′ \\ −α′ & 0 \end{matrix} ] , A_{j(r),j(r)}^{[r]+2} = [ \begin{matrix}0 & β′ \\ −β & 0 \end{matrix} ] .

From the relation (1.3) we have

∑_{r=0}^{K−1} ∥A_{i(r),j(r)}^{[r]}∥^{2} ≤ S(A)^{2}

So the Lemma follows by combining the above two inequalities.

In the sequel we also use i and j as indices for block-rows and blockcolumns and not as functions of r. In this context Φ_{ij} will be associated with A_{ij}^{[r]} where r is a unique number in {0, 1, …, K − 1} such that A_{ij}^{[r]} is the pivot submatrix of A^{[r]}. For example, the relation (2.13) implies that for the row-cyclic Paardekooper method holds

∑_{i=1}^{k−1} ∑_{j=i+1}^{k} Φ_{ij} ≤ [1.0062\frac{S(A)}{δ}]^{2} .  (2.14)

3 Quadratic Convergence

In this section we prove a sharp quadratic convergence bound for the rowcyclic Paardekooper method. In this section we heavily use the notations defined at the end of section 1. So the reader is advised to review the notations at this point.

Lemma 3.1 Suppose the row-cyclic Paardekooper method is applied to the skew-symmetric matrix A. Then after the completion of the first cycle holds

S(A^{[K]}) ≤ \frac{S(A)}{1 − max_{1≤i≤k−2}\{∑_{l=i+2}^{k} Φ_{il}\}}\sqrt{∑_{i=1}^{k−2} ∑_{l=i+2}^{k} Φ_{il}}.  (3.1)

Proof. First we estimate ∥A_{1j}^{[k−1]}∥^{2} for 2 ≤ j ≤ k − 1. Note that

| a_{pq}^{(r)} |≤ ω_{pq}, (p, q) ∈ Block(1, j), 4j − 4 ≤ r ≤ 2n − 4

where (see Figure 2 below)

ω_{1,2j−1} = ∑_{s=2j+1}^{n} |a_{2j−1,s}^{(1,s)}|| sin ϕ_{1s}| ω_{2,2j−1} = ∑_{s=2j+1}^{n} |a_{2j−1,s}^{(2,s)}|| sin ϕ_{2s}| ω_{1,2j} = ∑_{s=2j+1}^{n} |a_{2j,s}^{(1,s)}|| sin ϕ_{1s}| ω_{2,2j} = ∑_{s=2j+1}^{n} |a_{2j,s}^{(2,s)}|| sin ϕ_{2s}| .  (3.2)

If we set

Ω_{1j} = ∑_{(p,q)∈Block(1,j)} ω_{pq}^{2},  (3.3)

then

∥A_{1j}^{(r)}∥^{2} ≤ Ω_{1j}, 4j − 4 ≤ r ≤ 2n − 4, 2 ≤ j ≤ k − 1.  (3.4)

Using the Cauchy-Schwarz inequality to bound the right-hand sides of (3.2) and then using (3.3) we obtain for 2 ≤ j ≤ k − 1

Ω_{1j} ≤ ∑_{s=2j+1}^{n} |a_{2j−1,s}^{(1,s)}|^{2} · ∑_{s=2j+1}^{n} sin^{2} ϕ_{1s} + ∑_{s=2j+1}^{n} |a_{2j−1,s}^{(2,s)}|^{2} · ∑_{s=2j+1}^{n} sin^{2} ϕ_{2s}

+ ∑_{s=2j+1}^{n} |a_{2j,s}^{(1,s)}|^{2} · ∑_{s=2j+1}^{n} sin^{2} ϕ_{1s} + ∑_{s=2j+1}^{n} |a_{2j,s}^{(2,s)}|^{2} · ∑_{s=2j+1}^{n} sin^{2} ϕ_{2s}

≤ max_{1≤t≤2} { ∑_{s=2j+1}^{n} (|a_{2j−1,s}^{(t,s)}|^{2} + |a_{2j,s}^{(t,s)}|^{2})} · ∑_{l=j+1}^{k} Φ_{1l}.  (3.5)

To bound the term in braces in (3.5) we consider first the case n = 2k. From the order of rotations in ℑ_{1l} we observe (see Figure 2 below)

\begin{matrix}a_{r,2l}^{(2,2l)} & = & a_{r,2l}^{(1,2l−1)} & = & a_{r,2l}^{[1,l]} \\ a_{r,2l−1}^{(2,2l−1)} & = & a_{r,2l−1}^{(1,2l)} & = & a_{r,2l−1}^{[1,l]+2} \end{matrix} \} , r ∈ \{2j − 1, 2j\}

So we have

\begin{aligned}∑_{s=2j+1}^{n} (| a_{2j−1,s}^{(1,s)} |^{2} + | a_{2j,s}^{(1,s)} |^{2}) & = & ∑_{l=j+1}^{k} ; (∥A_{jl}^{[1,l]}e_{1}∥^{2} + ∥A_{jl}^{[1,l]+2}e_{2}∥^{2}), \\ ∑_{s=2j+1}^{n} (| a_{2j−1,s}^{(2,s)} |^{2} + | a_{2j,s}^{(2,s)} |^{2}) & = & ∑_{l=j+1}^{k} ; (∥A_{jl}^{[1,l]+2}e_{1}∥^{2} + ∥A_{jl}^{[1,l]}e_{2}∥^{2}), \end{aligned} (3.6)

where e_{1} = [1, 0]^{T}, e_{2} = [0, 1]^{T}.

[figure: 0; 0; 1; 2; 1; 1; 2; 1; 2j-1; 2j; j; ·; ·; ·; ·; ·; ·; ·; ·; ·; ·; ·; ·; ·; ·; ·; ·; x; x; x; x; x; x; x; x; A_{11}; A_{1j}; A_{1l}; A_{jj}; A_{jl}; pivot elements; affected elements; 2l-1; 2l; l; 0; 0; 2j-1; 2j; j]
Figure 2: j = 3, l = 6

To bound the elements of A_{jl}^{[1,l]+2} in terms of elements of A_{jl}^{[1,l]} and A_{1j}^{[1,l]} we use the obvious inequalities (see Figure 2)

|a_{2j−1,2l−1}^{[1,l]+2}| ≤ |a_{2j−1,2l−1}^{[1,l]}| + |a_{1,2j−1}^{[1,l]}|| sin ϕ_{1,2l−1}|
|a_{2j,2l−1}^{[1,l]+2}| ≤ |a_{2j,2l−1}^{[1,l]}| + |a_{1,2j}^{[1,l]}|| sin ϕ_{1,2l−1}|
|a_{2j−1,2l}^{[1,l]+2}| ≤ |a_{2j−1,2l}^{[1,l]}| + |a_{2,2j−1}^{[1,l]}|| sin ϕ_{2,2l}|
|a_{2j,2l}^{[1,l]+2}| ≤ |a_{2j,2l}^{[1,l]}| + |a_{2,2j}^{[1,l]}|| sin ϕ_{2,2l}|

Making use of the auxiliary inequality

(a + b)^{2} ≤ (1 + γ) a^{2} + (1 + \frac{1}{γ}) b^{2}, γ > 0  (3.7)

we obtain

∥A_{jl}^{[1,l]+2}e_{2}∥^{2} ≤ (1 + γ)∥A_{jl}^{[1,l]}e_{2}∥^{2} + (1 + \frac{1}{γ})∥e_{2}^{T} A_{1j}^{[1,l]}∥^{2} sin^{2} ϕ_{2,2l}
∥A_{jl}^{[1,l]+2}e_{1}∥^{2} ≤ (1 + γ)∥A_{jl}^{[1,l]}e_{1}∥^{2} + (1 + \frac{1}{γ})∥e_{1}^{T} A_{1j}^{[1,l]}∥^{2} sin^{2} ϕ_{1,2l−1}.

Thus

∥A_{jl}^{[1,l]}e_{1}∥^{2} + ∥A_{jl}^{[1,l]+2}e_{2}∥^{2} ≤ (1 + γ)∥A_{jl}^{[1,l]}∥^{2} + (1 + \frac{1}{γ})∥A_{1j}^{[1,l]}∥^{2} sin^{2} ϕ_{2,2l}

∥A_{jl}^{[1,l]+2}e_{1}∥^{2} + ∥A_{jl}^{[1,l]}e_{2}∥^{2} ≤ (1 + γ)∥A_{jl}^{[1,l]}∥^{2} + (1 + \frac{1}{γ})∥A_{1j}^{[1,l]}∥^{2} sin^{2} ϕ_{1,2l−1},

and the relations (3.5) and (3.6) yield

Ω_{1j} ≤[(1 + γ) ∑_{l=j+1}^{k} ∥A_{jl}^{[1,l]}∥^{2} + (1 + \frac{1}{γ}) max_{j+1≤l≤k} ∥A_{1j}^{[1,l]}∥^{2} ∑_{l=j+1}^{k} Φ_{1l}] · ∑_{l=j+1}^{k} Φ_{1l}.

(3.8)

One can verify that the relation (3.8) is still valid when n is odd.

By the relation (3.4) we have

∥A_{1j}^{[1,l]}∥^{2} ≤ Ω_{1j}, j + 1 ≤ l ≤ k.

So, the relation (3.8) implies

Ω_{1j} [1 − (1 + \frac{1}{γ})( ∑_{l=j+1}^{k} Φ_{1l})^{2}] ≤ (1 + γ) ∑_{l=j+1}^{k} ∥A_{jl}^{[1,l]}∥^{2} · ∑_{l=j+1}^{k} Φ_{1l}

which yields

Ω_{1j} ≤ \frac{1 + γ}{1 − (1 + \frac{1}{γ}) (∑_{l=j+1}^{k} Φ_{1l})^{2}} ∑_{l=j+1}^{k} ∥A_{jl}^{[1,l]}∥^{2} · ∑_{l=j+1}^{k} Φ_{1l}.

This inequality holds for all 2 ≤ j ≤ k − 1 and γ > 0. Since A_{1k}^{(2n−4)} = O we have

∑_{j=2}^{k} ∥A_{1j}^{(2n−4)}∥^{2} = ∑_{j=2}^{k−1} ∥A_{1j}^{(2n−4)}∥^{2} ≤ ∑_{j=2}^{k−1} Ω_{1j}

≤ \frac{1 + γ}{1 − (1 + \frac{1}{γ}) (∑_{l=3}^{k} Φ_{1l})^{2}} ∑_{j=2}^{k−1} ∑_{l=j+1}^{k} ∥A_{jl}^{[1,l]}∥^{2} · ∑_{l=3}^{k} Φ_{1l}.  (3.9)

To bound the double sum in the relation (3.9) we use invariance of the Euclidean norm under orthogonal transformations. Note that

∑_{j=1}^{l−1} ∥A_{jl}^{[1,l]}∥^{2} = ∑_{j=1}^{l−1} ∥A_{jl}^{[1,l−1]}∥^{2} = … = ∑_{j=1}^{l−1} ∥A_{jl}∥^{2}.

Hence

∑_{j=2}^{k−1} ∑_{l=j+1}^{k} ∥A_{jl}^{[1,l]}∥^{2} = ∑_{l=3}^{k} ∑_{j=2}^{l−1} ∥A_{jl}^{[1,l]}∥^{2} ≤ ∑_{l=3}^{k} ∑_{j=1}^{l−1} ∥A_{jl}^{[1,l]}∥^{2}

= ∑_{l=3}^{k} ∑_{j=1}^{l−1} ∥A_{jl}∥^{2} ≤ S^{2}(A)

and the relation (3.9) takes form

∑_{j=2}^{k} ∥A_{1j}^{[2,3]}∥^{2} = ∑_{j=2}^{k} ∥A_{1j}^{(2n−4)}∥^{2} ≤ \frac{(1 + γ)S^{2}(A)}{1 − (1 + \frac{1}{γ}) (∑_{l=3}^{k} Φ_{1l})^{2}} ∑_{l=3}^{k} Φ_{1l}, γ > 0.

Applying the same analysis to the i’th block-row we obtain for 1 ≤ i ≤ k −2

∑_{j=i+1}^{k} ∥A_{ij}^{[i+1,i+2]}∥^{2} ≤ \frac{1 + γ}{1 − (1 + \frac{1}{γ}) (∑_{l=i+2}^{k} Φ_{il})^{2}} S^{2}(A^{[i,i+1]}) ∑_{l=i+2}^{k} Φ_{il}.

In all inequalities we have used the same γ > 0. By (1.3) S^{2}(A^{[i,i+1]}) ≤ S^{2}(A), 1 ≤ i ≤ k − 2. Note that the sum ∑_{j=i+1}^{k} ∥A_{ij}^{[i+1,i+2]}∥^{2} is invariant till the end of the current cycle and because of the choice of pivot strategy A_{k−1,k}^{[K]} = O. Hence. for γ > 0 holds

S^{2}(A^{[K]}) = ∑_{i=1}^{k−2} ∑_{j=i+1}^{k} ∥A_{ij}^{[i+1,i+2]}∥^{2} ≤ \frac{(1 + γ)S^{2}(A) ∑_{i=1}^{k−2} ∑_{l=i+2}^{k} Φ_{il}}{1 − (1 + \frac{1}{γ}) max_{1≤i≤k−2} (∑_{l=i+2}^{k} Φ_{il})^{2}} .

(3.10)

To minimize the right-hand side of (3.10) with respect to γ we consider the function

f(γ) = \frac{1 + γ}{1 − (1 + \frac{1}{γ} ) · x}, γ > 0, x ≥ 0.

The function f is minimum at γ_{0} = \sqrt{x}/(1 − \sqrt{x}) and f(γ_{0}) = 1/(1 − \sqrt{x})^{2}. With x = max_{1≤i≤k−2} (∑_{l=i+2}^{k} Φ_{il})^{2} we obtain (3.1).

Theorem 3.2 Let a skew-symmetric matrix A sutisfy the condition (2.2). If the row-cyclic Paardekooper method is applied to A, then after the completion of the first cycle we have

S(A^{[K]}) ≤ 1.35\frac{S^{2}(A)}{δ}.  (3.11)

The inequality (3.11) holds with any cyclic strategy which is equivalent to the row-cyclic one.

Proof. Since max_{1≤i≤k−2}{∑_{l=i+2}^{k} Φ_{il}} ≤ ∑_{i=1}^{k−2} ∑_{l=i+2}^{k} Φ_{il}, Lemma 3.1 implies

S(A^{[K]}) ≤ \frac{S(A)}{1 − ∑_{i=1}^{k−2} ∑_{l=i+2}^{k} Φ_{il}} \sqrt{∑_{i=1}^{k−2} ∑_{l=i+2}^{k} Φ_{il}}.  (3.12)

From the relation (2.14) and the assumption (2.2) we obtain

∑_{i=1}^{k−2} ∑_{l=i+2}^{k} Φ_{il} ≤ 1.0062^{2} \frac{S^{2}(A)}{δ^{2}} < 0.5031^{2}.  (3.13)

From the relations (3.12) and (3.13) we obtain

S(A^{[K]}) ≤ \frac{S(A)}{1 − 0.5031^{2}} · 1.0062 · \frac{S(A)}{δ} ≤ 1.35\frac{S^{2}(A)}{δ}

which completes the proof for the row-cyclic strategy. Since the same matrix A^{[K]} is obtained for any cycle strategy which is equivalent to the row-cyclic one (see [3, 4]) the theorem is proved.

With another measure

∥A∥_{OF F} = \sqrt{2}S(A) = [∑_{i=1}^{k} ∑_{j≠i^{j=1}}^{k} ∥A_{ij}∥^{2}]^{\frac{1}{2}} ,

Theorem 3.2 can be written in the following simple form

∥A∥_{OF F} < \frac{\sqrt{2}}{2}δ ⟹ ∥A^{[K]}∥_{OF F} ≤ \frac{∥A∥_{OF F}^{2}}{δ}

This sharp quadratic convergence bound of the row-cyclic Paardekooper method should be compared to that of a general cyclic Paardekooper method in [7], which says

∥A∥_{OF F} < \frac{\sqrt{2}}{2}δ ⟹ ∥A^{[K]}∥_{OF F} ≤ 0.5668(n − 1)\frac{∥A∥_{OF F}^{2}}{δ}

From Lemma 3.1 we see that Φ_{i,i+1}, 1 ≤ i ≤ k − 1 do not enter into the bound for S(A^{[K]}). This is due to the special cyclic, that is row-cyclic, strategy. Using that fact we can relax the assumptions of Theorem 3.2. Let us remember that each diagonal block A_{jj} is affiliated with a quantity ν_{j} where ±ıν_{j} is the eigenvalue pair of A. In order to allow double eigenvalues or clusters of two simple eigenvalues we propose another measure δ̄ of place of δ. It is computed via the following algorithm.

Algorithm 3.3 Input data are real numbers ν_{j}, 1 ≤ j ≤ k. Output is δ̄.

1. Sort ν_{1}, …, ν_{k} in the nonincreasing order

2. Compute μ:

2.1 μ = ν_{1}

2.2 for i = 2, …, k − 1 do γ = max{ν_{i−1} − ν_{i} , ν_{i} − ν_{i+1}} if γ < μ then μ = γ endfor

3. δ̄ = μ/3.

Note that δ̄ > 0 implies at most double (if ν_{1} = 0 at most quadruple) eigenvalues. It is obvious (see Figure 3) that δ̄ can be significantly larger than δ. And we always have δ̄ ≥ δ.

[figure: ^{3δ̄}; ν_{7}; ν_{6}; ν_{5}; ν_{4} = ν_{3}; ν_{2}; ν_{1}; 0]

Figure 3: The quantity δ̄ If δ̄ > 0, we define

Cl_{i} = {ν_{j}; |ν_{j} − ν_{i}| < 3δ̄, 1 ≤ j ≤ k}, 1 ≤ i ≤ k.

Note that ν_{i} ∈ Cl_{i}, i = 1, …, k. Hence the sets Cl_{i}, i = 1, …, k are nonempty sets, and they are either identical or disjoint and each of them can contain at most two (identical or distinct) elements. We call Cl_{i} the cluster of the eigenvalue ν_{i}, i = 1, …, k. Using δ̄ we propose the first asymptotic assumption

A1: S(A) < \frac{δ̄}{2}, δ̄ > 0.

Note that the assumption A1 implies δ̄ > 0. Hence each Cl_{i}, i = 1, …, k, has at most two elements including ν_{i}. As before we say that the diagonal block A_{ii}, i = 1, …, k′, is affiliated with ν_{j} if ||a_{2i−1,2i}| − ν_{j}| < δ̄/2. Due to Hoffman–Wielandt theorem, we may assume that A_{ii} is affiliated with ν_{i}, i = 1, …, k′. If ν_{j} ∈ Cl_{i}, we say that A_{jj} is associated with A_{ii}. Since Cl_{i} has at most two elements, at most two diagonal blocks are associated. Now we propose the second asymptotic assumption.

A2: The associated diagonal blocks occupy successive positions on the blockdiagonal.

We assume that the given skew-symmetric A satisfies both assumptions A1 and A2, and we use the row-cyclic Paardekooper method. The proof of Lemma 2.1 shows that the affiliation of diagonal blocks can be changed only if two diagonal blocks are associated. This means that the association of diagonal blocks is preserved. Hence all A^{[r]}, r ≥ 0, satisfies the assumption A2. A simple inspection shows that the angle estimates of Lemma 2.2 with δ̄ in place of δ hold provided that A_{ii} and A_{jj} are not associated. So for Lemma 2.3 to hold we have to replace ∑_{r=0}^{K−1} Φ_{i(r),j(r)} by ∑′Φ_{i(r),j(r)} where ∑′ denotes the sum over r ∈ {0, …, K − 1} for which A_{i(r),i(r)}^{[r]} and A_{j(r),j(r)}^{[r]} are not associated. That is, instead of (2.13) we have

∑′Φ_{i(r),j(r)} ≤ [1.0062\frac{S(A)}{δ̄}]^{2}  (3.14)

Since the number of blocks which are associated is at most two and the assumption A2 is satisfied by all A^{[r]}, r ≥ 0, those associated diagonal blocks are next to each other on the block-diagonal. Hence it follows that

∑_{i=1}^{k−2} ∑_{l=i+2}^{k} Φ_{il} ≤ ∑′Φ_{i(r),j(r)}  (3.15)

Now by following the proof of Theorem 3.2 and using the relations (3.14) and (3.15) and the assumption A1 we obtain the following Theorem which is an improvement of Theorem 3.2.

Theorem 3.4 Let the assumptions A1 and A2 hold. Then after the completion of the first cycle of the row-cyclic Paardekooper method we have

S(A^{[K]}) ≤ 1.35\frac{S^{2}(A)}{δ̄}.

This inequality holds for Paardekooper method under any cyclic strategy which is equivalent to the row-cyclic one.

Numerical Examples

In the following we will consider some numerical examples. All the computations were done on VAX6000 Model 500 using quadruple precision arithmetic with machine precision eps ≃ 9.63 x 10^{−35}.

Example 3.5 We have generated matrix A of order 40 from the Murnaghan form whose nontrivial entries are 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10 by a series of 2n^{2} similarity transformations with plane rotations. For each matrix generation we have used the row-cyclic ordering and rotation angles are chosen randomly from the interval [−\frac{π}{4} , \frac{π}{4}]. In this example δ̄ = 1/3 and δ = 0 and it has been checked that A satisfies the assumptions A1 and A2.

By the same procedure we have generated matrix Ã of order 40 with eigenvalues 1, 2, 1, 2, 3, 4, 3, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10. Here δ̄ = 1/3, δ = 0 and it has been checked that Ã satisfies the assumptions A1, but not the assumption A2.

Applying the row-cyclic Paardekooper method to A and Ã we have computed

\frac{S(A^{[tK]})}{δ̄}, s_{t} = \frac{S(A^{[tK]})/δ̄}{1.35[S(A^{[(t−1)K]})/δ̄]^{2}} , \frac{S(Ã^{[tK]})}{δ̄}, s̃_{t} = \frac{S(Ã^{[tK]})/δ̄}{1.35[S(Ã^{[(t−1)K]})/δ̄]^{2}} .

Here t counts cycles of the processes. By Theorem 3.4 s_{t} cannot be larger than 1. We have obtained the following results:

cycle  S(A)/δ̄  s_{t}  S(Ã)/δ̄  s̃_{t}
t=0 1.814 · 10^{−1}  1.811 · 10^{−1}
t=1 1.941 · 10^{−4} 4.368 · 10^{−3} 5.619 · 10^{−3} 1.269 · 10^{−1}
t=2 4.229 · 10^{−10} 8.317 · 10^{−3} 4.645 · 10^{−5} 1.090 · 10^{0}
t=3 3.934 · 10^{−21} 1.630 · 10^{−2} 1.871 · 10^{−9} 6.426 · 10^{−1}
t=4 3.428 · 10^{−43} 1.640 · 10^{−2} 1.496 · 10^{−17} 3.164 · 10^{0}
t=5 3.460 · 10^{−87} 2.181 · 10^{−2} 2.039 · 10^{−18} 6.751 · 10^{15}
t=6 1.740 · 10^{−175} 1.077 · 10^{−2} 1.672 · 10^{−23} 2.979 · 10^{12}

We note the quadratic reduction of S(A)/δ̄ after each cycle. By looking at s_{t} column we see that the quadratic convergence bound in Theorem 3.4 is only about 10^{2} times larger than the actual constant. This size of overestimate seems to be independant of the size of S(A)/δ̄. Similar results were obtained for matrices with simple eigenvalues satisying the assumption A1. We also note that S(Ã)/δ̄ fails to converge quadratically (for t = 2, 4, 5, 6, s̃_{t} is larger than 1) because Ã does not satisfy the assumption A2.

In the above example we have considered double eigenvalues. In the following example we consider clusters of two eigenvalues.

Example 3.6 By the same procedure as in Example 3.5 we have generated matrices A and Ã of order 41 from the Murnaghan forms whose nontrivial entries are 1, 1 + 10^{−15}, 3, 3 + 10^{−15}, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12 and 1, 3, 1+10^{−15}, 3+10^{−15}, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, respectively. For the both matrices we have δ = 0 and δ̄ = 1/3. It has been checked that the initial matrix A satisfies the assumptions A1 and A2, but for Ã the assumption A2 is not satisfied, even though A1 is satisfied. Applying the row-cyclic Paardekooper method to A and Ã we have obtained the following results:

cycle  S(A)/δ̄  s_{t}  S(Ã)/δ̄  s̃_{t}
t=0 2.434 · 10^{−1}  2.433 · 10^{−1}
t=1 2.235 · 10^{−4} 2.794 · 10^{−3} 4.882 · 10^{−3} 6.110 · 10^{−2}
t=2 5.173 · 10^{−10} 7.668 · 10^{−3} 1.051 · 10^{−5} 3.268 · 10^{−1}
t=3 2.688 · 10^{−21} 7.440 · 10^{−3} 6.828 · 10^{−9} 4.576 · 10^{1}
t=4 1.714 · 10^{−43} 1.758 · 10^{−2} 7.108 · 10^{−10} 1.129 · 10^{7}
t=5 3.650 · 10^{−88} 9.202 · 10^{−3} 5.096 · 10^{−17} 7.472 · 10^{1}
t=6 6.530 · 10^{−177} 3.630 · 10^{−2} 2.050 · 10^{−40} 5.846 · 10^{−8}

We note the quadratic reduction of S(A)/δ̄ after each cycle. We note a delay of the quadratic convergence of S(Ã)/δ̄ (look at s̃_{t} column when t = 3, 4, 5). It is caused because Ã does not satisfy the assumption A2. Ultimately the quadratic reduction again prevails since 1, 3, 1+10^{−15}, 3+10^{−15} cannot be treated as clustered eigenvalues when S(Ã) < \frac{1}{6} 10^{−15}.

Our next example shows that in the case of eigenvalues with multiplicities larger than two (or clusters with more than two eigenvalues) we cannot expect the quadratic convergence. This is in sharp contrast with the case of symmetric Jacobi method. For the symmetric Jacobi method, if we let 3δ be the minimal distance between two distinct eigenvalues then the row cyclic Jacobi method yields ultimate quadratic reduction of S(A)/δ (see [6]) .

Example 3.7 We have generated a matrix of order 16 from the Murnaghan form whose nontrivial entries are 1, 1, 1, 2, 2, 2, 3, 4 by the same procedure as in the previous examples. Here δ′ can be defined as 3δ′ = max_{ν_{i}≠ν_{j}} |ν_{i} − ν_{j}| and thus δ′ = 1/3. It has been checked that the initial matrix A satisfies the assumptions A1 (with δ′, not δ̄) and A2. By the same procedure we have generated matrix Ã of order 16 with eigenvalues 1 − 2 · 10^{−15}, 1 − 10^{−15}, 1, 2, 2 + 10^{−15}, 2 + 2 · 10^{−15}, 4, 4. Since eigenvalues are clustered near 1, 2 and 4, we define 3δ̃ to be the smallest distance between two clusters. Thus δ̃ = 1/3. It has been checked that the initial matrix Ã satisfies the assumptions A1 (with δ̃, not δ̄) and A2. Applying the rwo-cyclic Paardekooper method to A and Ã we have obtained the following results:

cycle S(A)/δ′  s_{t}  S(Ã)/δ̃  s̃_{t}
t=0 2.431 · 10^{−1}  2.606 · 10^{−2}
t=1 5.486 · 10^{−3} 6.877 · 10^{−2} 4.023 · 10^{−4} 4.387 · 10^{−1}
t=2 1.136 · 10^{−4} 2.797 · 10^{0} 1.097 · 10^{−7} 5.020 · 10^{−1}
t=3 5.089 · 10^{−10} 2.920 · 10^{−2} 1.690 · 10^{−9} 1.040 · 10^{5}
t=4 2.654 · 10^{−19} 7.590 · 10^{−1} 1.866 · 10^{−11} 4.842 · 10^{6}
t=5 2.651 · 10^{−20} 2.788 · 10^{17} 1.160 · 10^{−18} 2.467 · 10^{3}
t=6 3.179 · 10^{−22} 3.350 · 10^{17} 5.729 · 10^{−49} 3.153 · 10^{−13}

We note that S(A)/δ fails to converge quadratically (look at s_{t} column when t = 5, 6). It is due to the fact that A has eigenvalues whose multiplicities are larger than 2. We also note a delay of the quadratic reduction of S(Ã)/δ̃ (look at s̃_{t} column when t = 3, 4, 5). It is because Ã has clusters of more than two simple eigenvalues. Ultimately the quadratic reduction of S(Ã)/δ̃ again prevails since 1 − 2 · 10^{−15}, 1 − 10^{−15}, 1, 2, 2 + 10^{−15}, 2 + 2 · 10^{−15} cannot be treated as clustered eigenvalues when S(Ã) < \frac{1}{6} 10^{−15}.

The last example shows that δ̄ defined by Algorithm 3.3 is the best possible one under which the quadratic convergence of S(A)/δ̄ is guaranteed.

References

[1] J. Demmel and K. Veselić, Jacobi’s Method is More Accurate than QR”, Courant Institute of Mathematical Sciences, Department of Computer Science, New York University 251 Mercer Street, New York, Y. Y. 10012, October, 1989.

[2] K. K. Gupta, On a Combined Strum Sequence and Inverse Iteration Technique for Eigenproblem of Spinning Structures, International Journal for Numerical Methods in Engineering 7: 509–518 (1973).

[3] E. R. Hansen, On Jacobi Methods and Block Jacobi Methods for Computing Matrix Eigenvalues, Ph. D. thesis, Stanford University 1960.

[4] E. R. Hansen, On Cyclic Jacobi Methods, SIAM J. Appl. Math. 11: 449–459 (1963).

[5] V. Hari, Convergence of the Paardekooper Method for Diagonalization of Skew-symmetric matrices, M.S.Thesis, University of Zagreb, 1980 (Croatian language).

[6] V. Hari, On Sharp Quadratic Convergence Bounds for the Serial Jacobi Methods, Numer. Math. 60: 375–406 (1991).

[7] V. Hari, On Quadratic Convergence of the Paardekooper Method, Glasnik Matematički 17 (37): 183–195 (1982).

[8] A. J. Hoffman and H.W. Wielandt, The Variation of the Spectrum of a Normal Matrix, Duke Math. J. 20: 37-39 (1953).

[9] F. D. Murnaghan and A. Winter, A Canonical Form for Real Matrices under Orthogonal Transformations, USA Nat. Acad. Sci. 17: 417-420 (1931).

[10] M. H. C. Paardekooper, An Eigenvalue Algorithm for Skew-symmetric Matrices, Numer.Math. 17: 189-202 (1971).

[11] N.H.Rhee and V.Hari, On the Global and Cubic Convergence of a Quasicyclic Jacobi Method, Technical Report 91-27, AHPCRC, University of Minnesota.

[12] R. A. Rosanoff et al. , Numerical Conditions of Stiffness Matrix Formulations for Frame Structures, Proceedings of the 2nd Conference on Matrix Methods in Structural Mechanics, WPAFB Dayton, Ohio, 1968.

[13] J. H. Wilkinson, Note on the Quadratic Convergence of the Cyclic Jacobi Process, Numerische Mathematik 4: 296–300 (1962).