ON THE CUBIC CONVERGENCE OF THE PAARDEKOOPER METHOD

NOAH H. RHEE and VJERAN HARI

Department of Mathematics Department of Mathematics University of Missouri-Kansas City University of Zagreb Kansas City, MO 64110-2499 Bijenička cesta 30, 41001 Zagreb USA Croatia e-mail: rhee@vax1.umkc.edu e-mail: hari@math.hr

Abstract.

We show a simple way how asymptotic convergence results can be conveyed from a simple Jacobi method to a block Jacobi method. Our pilot methods are the well known symmetric Jacobi method and the Paardekooper method for reducing a skew-symmetric matrix to the real Schur form. We show resemblance in the quadratic and cubic convergence estimates, but also discrepances in the asymptotic assumptions. By numerical tests we confirm that our asymptotic assumptions for the Paardekooper method are most general.

AMS Subject Classifications: 65F15.

Keywords: Skew-symmetric matrix, Paardekooper method, cubic convergence.

1. Introduction

In this paper we derive a sharp cubic convergence bound of the Paardekooper method for skewsymmetric matrices under a suitable quasi-cyclic strategy. Similar as in [12] we deal with the row-cyclic representative from the Mascarenhas class [9] of cubically convergent quasi-cyclic pivot strategies. Even though the Paardekooper method operates on 2-by-2 blocks, the analysis is much the same as that of the symmetric Jacobi method [12]. We tried to use as much as possible the results from [12]. Hence the obtained asymptotic estimates are very similar to those in the symmetric case. However, there is a difference in the asymptotic conditions. This is a consequence of the special property of a skew-symmetric matrix which is almost in its real Schur (i.e. Murnaghan) form [7].

The skew-symmetric eigenvalue problems appear in different areas of applied mathematics [1, 13]. The Paardekooper method is a reliable diagonalization method, suitable for parallel and vector computations. It can also be efficiently implemented on standard machines for which simple 4-by-4 matrix transformation instructions are available (for example, IIT coprocessor supported software).

The paper is organized as follows. In Section 2 we discuss a Paardekooper step and in Section 3 we define the quasi-cyclic Paardekooper method. In Section 4 we improve the quadratic convergence estimate from [6] under the row-cyclic strategy. In section 5 we establish the cubic convergence of the quasi-cyclic method and present numerical tests.

2. Paardekooper Step

Let A = (a_{ij}) be a skew-symmetric matrix of order n = 2k. If n is odd we can append to A one zero row and one zero column. By A = (A_{ij}) , 1 ≤ i, j ≤ k, we denote the Paardekooper partition [11] of A, where each submatrix A_{ij} is 2-by-2. The Paardekooper method is a Jacobi-type method which repeatedly annihilates the off-diagonal blocks A_{ij}, thus reducing the measure

S (A) = \sqrt{∑_{i<j} ∥A_{ij}∥_{F}^{2}} ,

and bringing the matrix closer to the Murnaghan form [10]. Here ∥ · ∥_{F} stands for the Euclidean (Frobenius) matrix norm. One step of the method consists of the orthogonal similarity transformation

A′ = J_{pq}^{T}AJ_{pq}, p < q,

which makes A′_{pq} zero. This ensures that

S^{2} (A′) = S^{2} (A) − ∥A_{pq}∥_{F}^{2} (1)

holds. The matrix J_{pq} is referred to as the Jacobi annihilator of A_{pq}, A_{pq} is the pivot submatrix of A and (p, q) is the pivot pair. Let R (i, j; ϕ) stand for the rotation matrix in the (i, j) plane with angle ϕ, whose (i, i), (i, j), (j, i) and (j, j) elements are cos ϕ, sin ϕ, − sin ϕ and cos ϕ, respectively. The Jacobi annihilator J_{pq} is built of four plane rotations:

J_{pq} = R (2p − 1, 2q − 1; ϕ_{1}) R (2p, 2q; ϕ_{2}) R (2p − 1, 2q; ϕ_{3}) R (2p, 2q − 1; ϕ_{4}) . (2)

Using simplified notation

A_{pp} = [ \begin{matrix}0 & α \\ −α & 0 \end{matrix} ] , A_{pq} = [ \begin{matrix}x & v \\ w & y \end{matrix} ] , A_{qq} = [ \begin{matrix}0 & β \\ −β & 0 \end{matrix} ] ,

the angle formulas (cf. [11, 3, 5]) for ϕ_{1} ∈ (−π/4, π/4] and ϕ_{2} ∈ (−π/2, π/2] are

tan 2ϕ_{1} = \frac{2 (αw − βv)}{α^{2} − β^{2} + v^{2} − w^{2}} , tan ϕ_{2} = \{ \begin{matrix}−\frac{c_{1}v−s_{1}β}{c_{1}α+s_{1}w} , & c_{1}α + s_{1}w ≠ 0 \\ −\frac{c_{1}w−s_{1}α}{c_{1}β+s_{1}v}, & c_{1}β + s_{1}v ≠ 0 \end{matrix} . (3)

Here c_{1} = cos ϕ_{1} and s_{1} = sin ϕ_{1}. The first two angles have been chosen to annihilate v and w. Note that x and y have not changed. Suppose that A_{pp}, A_{pq} and A_{qq} are transformed by the first two rotations into

Ã_{pp} = [ \begin{matrix}0 & α̃ \\ −α̃ & 0 \end{matrix} ] , Ã_{pq} = [ \begin{matrix}x & 0 \\ 0 & y \end{matrix} ] , Ã_{qq} = [ \begin{matrix}0 & β̃ \\ −β̃ & 0 \end{matrix} ]

The angles ϕ_{3} ∈ (−π/4, π/4] and ϕ_{4} ∈ (−π/2, π/2] are computed from the requirement that x and y become zero, which yields

tan 2ϕ_{3} = \frac{2 (α̃y + β̃x)}{α^{2} − β̃^{2} + x^{2} − y^{2}} , tan ϕ_{4} = \{ \begin{matrix}−\frac{c_{3}x+s_{3}_{β}}{c_{3}_{α}+s_{3}y}, & c_{3}α̃ + s_{3}y ≠ 0 \\ \frac{c_{3}y−s_{3}_{α}}{c_{3}_{β}−s_{3}x}, & c_{3}β̃ − s_{3}x ≠ 0 \end{matrix} , (4)

where c_{3} = cos ϕ_{3} and s_{3} = sin ϕ_{3}. For a more detailed description of the Jacobi annihilator J_{pq} see [11, 5, 6].

Now consider the case when originally A is of order n = 2k − 1 and we have appended to it the zero row and column. Then in A_{pk} and A_{kk}, v = y = 0 and β = 0, respectively, and the formulas (3) and (4) imply

tan 2ϕ_{1} = \frac{2αw}{α^{2} − w^{2}} , tan ϕ_{2} = 0, tan 2ϕ_{3} = 0, tan ϕ_{4} = −\frac{x}{α̃}.

Hence

J_{pk} = R (2p − 1, 2k − 1; ϕ_{1}) R (2p, 2k − 1; ϕ_{4})

is made of just two instead of four plane rotations. This J_{pk} is slightly different from the corresponding Jacobi annihilator given in [11, 3, 5, 6].

Since R (2p, 2k; ϕ_{2}) and R (2p − 1, 2k; ϕ_{3}) are the identities, it is easy to see that the zero structure of the last row and column of A is preserved by the transformation with J_{pk}. This means that we can actually work with the original matrix A, but for the sake of analysis of the algorithm we further assume n = 2k.

3. The Quasi-Cyclic Paardekooper Method P_{Θ}

To completely describe a Paardekooper method we need to define a pivot strategy. We identify pivot strategy with function Π from N = {1, 2, ⋯} to P = {(p, q) : 1 ≤ p < q ≤ k}. It means that Π (r) = (p (r) , q (r)) denotes the pivot pair at step r. The Paardekooper method, defined by a pivot strategy Π, is denoted by P (Π) . If Π is a periodic function, we call it a periodic (pivot) strategy. Note that K = k (k − 1) /2 is the cardinallity of the set P. Let Π be a periodic strategy with period M. If M > K (M = K) and {Π (r) : r = 1, ⋯, M} = P, Π is called a quasi-cyclic (cyclic) strategy. Every sequence of M successive steps starting at position Π(1) is called quasi-sweep (sweep) of the method. If Π is a quasi-cyclic (cyclic) strategy, P (Π) is called a quasi-cyclic (cyclic) Paardekooper method.

Next we define the quasi-cyclic strategy we deal with in this paper (cf. [12]). First we introduce some notation. Let S be a subset of P. By Ω_{R} (S) we denote the sequence which represents the row-wise ordering of S. That is, the elements of S in Ω_{R} (S) are ordered to satisfy the following condition: (p_{1}, q_{1}) ∈ S precedes (p_{2}, q_{2}) ∈ S iff either p_{1} < p_{2} or p_{1} = p_{2} and q_{1} < q_{2}. In an obvious way one can define Ω_{C} (S), the column-wise ordering of S. If S_{1}, ⋯, S_{α} are subsets of P, by Ω_{R} (S_{1}) , ⋯, Ω_{R} (S_{α}) we mean the sequence which is obtained by the concatination of the sequences Ω_{R} (S_{1}) , ⋯, Ω_{R} (S_{α}).

Let k_{1}, k_{2}, k_{3} and k_{4} be positive integers satisfying k_{1} + k_{2} + k_{3} + k_{4} = k and let

P_{1} = {(p, q) : 1 ≤ p < q ≤ k_{1}}
P_{2} = {(p, q) : k_{1} + 1 ≤ p < q ≤ k_{1} + k_{2}}
P_{3} = {(p, q) : k_{1} + k_{2} + 1 ≤ p < q ≤ k_{1} + k_{2} + k_{3}}
P_{4} = {(p, q) : k_{1} + k_{2} + k_{3} + 1 ≤ p < q ≤ k}
P_{5} = {(p, q) : 1 ≤ p ≤ k_{1} , k_{1} + 1 ≤ q ≤ k_{1} + k_{2}}
P_{6} = {(p, q) : k_{1} + k_{2} + 1 ≤ p ≤ k_{1} + k_{2} + k_{3} , k_{1} + k_{2} + k_{3} + 1 ≤ q ≤ k}
P_{7} = {(p.q) : 1 ≤ p ≤ k_{1} + k_{2} , k_{1} + k_{2} + 1 ≤ q ≤ k}.

Then {P_{i}, 1 ≤ i ≤ 7} makes a partition of P. With this partition of P is closely related the following block partition of A :

A = [ \begin{matrix}A_{1} & A_{3} \\ −A_{3}^{T} & A_{2} \end{matrix} ] , A_{1} = [ \begin{matrix}B_{11} & B_{12} \\ −B_{12}^{T} & B_{22} \end{matrix} ] , A_{2} = [ \begin{matrix}B_{33} & B_{34} \\ −B_{34}^{T} & B_{44} \end{matrix} ] . (5)

Here B_{11}, B_{22}, B_{33}, and B_{44} are square matrices of order 2k_{1}, 2k_{2}, 2k_{3}, and 2k_{4}, respectively. Note that (p, q) ∈ P_{i}, 1 ≤ i ≤ 4 iff A_{pq} belongs to the upper triangular part of B_{ii}. Also, (p, q) ∈ P_{5} (P_{6}) iff A_{pq} belongs to B_{12} (B_{34}) and (p, q) ∈ P_{7} iff A_{pq} belongs to A_{3}. We set Ω_{i} = Ω_{R} (P_{i}) , 1 ≤ i ≤ 7.

Definition 1. Let

Θ = Ω_{1}, Ω_{2}, Ω_{3}, Ω_{4}, Ω_{5}, Ω_{6}, Ω_{1}, Ω_{2}, Ω_{3}, Ω_{4}, Ω_{7}

and let M be the number of pairs in the sequence Θ. Let Π_{Θ} be the quasi-cyclic strategy with period M generated by Θ. The corresponding quasi-cyclic Paardekooper method P (Π_{Θ}) is denoted by P_{Θ}. □

If one assumes Ω_{i} = Ω_{C} (P_{i}), 1 ≤ i ≤ 7 in Definition 1, then one obtains a method which is essentially the same as P_{Θ} (cf. [12, rel. (2.1)],[4]). In this paper we discuss the ultimate cubic reduction of S (A) per quasi-sweep of the method P_{Θ}.

4. The Row-Cyclic Paardekooper Method

To analyze the effect of Paardekooper annihilators J_{pq} when (p, q) assumes the values from the sequences Ω_{1}, Ω_{2}, Ω_{3}, and Ω_{4}, we need to study the row-cyclic Paardekooper method. The quadratic convergence of the row-cyclic Paardekooper method has been studied in our earlier paper [6]. Here we explain the techique how the estimates from the symmetric Jacobi method can be conveyed to the Paardekooper method. As a result of that we improve the result given in [6].

Let A^{[0]} = A be of order n = 2k and let A^{[r]} = (A_{ij}^{[r]}) be the matrix obtained after the rth step of the row-cyclic Paarderkooper method. Since all the steps we consider belong to one sweep we denote by ϕ_{ij} the angle used in the (i, j) plane. By A^{[p,q]} we denote the matrix obtained just before the transformation with J_{pq} takes place. We also define

Block (p, q) = {(2p − 1, 2q − 1) , (2p, 2q) , (2p − 1, 2q) , (2p, 2q − 1)}

and

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

Note that the transformation with J_{pq} couples together the blocks A_{ip} and A_{iq,} (A_{pi} and A_{qi}), 1 ≤ i ≤ n, i ≠ p, q. For these coupled blocks we have the following

Lemma 2. For 1 ≤ i ≤ k, i ≠ p, q holds

∥A′_{ip}∥_{F} ≤ ∥A_{ip}∥_{F} + ∥A_{iq}∥_{F} Φ_{pq} , ∥A′_{iq}∥_{F} ≤ ∥A_{iq}∥_{F} + ∥A_{ip}∥_{F} Φ_{pq}.

Proof. . Note that

[A′_{ip}, A′_{iq}] = [A_{ip}, A_{iq}] Ĵ, i ≠ p, q

holds, where Ĵ stands for the 4-by-4 restriction of J_{pq} to the subspace spanned by the coordinate vectors e_{2p−1}, e_{2p}, e_{2q−1} and e_{2q}. Using the relation (2) we obtain

Ĵ = [ \begin{matrix}J_{11} & J_{12} \\ J_{21} & J_{22} \end{matrix} ] = [ \begin{matrix}c_{1}c_{3} & −s_{1}s_{4} & s_{1}c_{4} & s_{3}c_{1} \\ −s_{2}s_{3} & c_{2}c_{4} & s_{4}c_{2} & s_{2}c_{3} \\ −s_{1}c_{3} & −s_{4}c_{1} & c_{1}c_{4} & −s_{1}s_{3} \\ −s_{3}c_{2} & −s_{2}c_{4} & −s_{2}s_{4} & c_{2}c_{3} \end{matrix} ] , \begin{matrix}c_{i} = cos ϕ_{i} \\ s_{i} = sin ϕ_{i} \end{matrix} , 1 ≤ i ≤ 4,

where each J_{ij} is a 2-by-2 submatrix of J. Then obviously, A′_{ip} = A_{ip}J_{11} + A_{iq}J_{21} and it follows

∥A′_{ip}∥_{F} ≤ ∥A_{ip}∥_{F} ∥J_{11}∥_{2} + ∥A_{iq}∥_{F} ∥J_{21}∥_{F} .

Here ∥·∥_{2} denotes the spectral matrix norm. Since Ĵ is orthogonal and the spectral norm of any principal submatrix is not greater than the norm of the original matrix, ∥J_{11}∥_{2} ≤ 1. An easy calculation shows that ∥J_{21}∥_{F} ≤ Φ_{pq}. The second assertion is obtained in a similar way. □

Because of the skew-symmetry of A′ and A, in the assertions of Lemma 2, we can replace (independently) A′_{ip} by A′_{pi} (A′_{iq} by A′_{qi} ) and A_{ip} by A_{pi} (A_{iq} by A_{qi}).

We note that a similar analysis was done by Hacon [2] for a quaternion Jacobi method for skew symmetric matrices. The following lemma improves the estimate in [6, Lemma 3.1] and thus decreases the bound 4/3 in [6, Theorem 3.4] to 1.

Lemma 3. Suppose the row-cyclic Paardekooper method is applied to A. Then after the completion of the first sweep we have

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

Proof. . Applying the analysis of Wilkinson from [14] to the first block row of A and using Lemma 2 we obtain

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

Squaring, using the Cauchy-Schwarz inequality and summing over 2 ≤ j ≤ k − 1 we obtain

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

Using the invariance of the Frobenius norm under orthogonal transformations (analysis is the same as in [14]) we can bound the double sum in the relation (6) by S^{2} (A). Since A^{[k−1]} = A^{[2,3]}, and A_{1k}^{[k−1]} is is the null matrix, the relation (6) implies

∑_{j=2}^{k} ∥A_{1j}^{[2,3]}∥_{F}^{2} ≤ S^{2} (A) ∑_{l=3}^{k} Φ_{1l}^{2}.

Applying the same analysis to the ith block-row of A^{[i,i+1]}, we obtain

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

By the relation (1) we have S^{2} (A^{[i,i+1]}) ≤ S^{2} (A) , 1 ≤ i ≤ k−2. Note that the sum ∑_{j=i+1}^{k} ∥A_{ij}^{[t]}∥_{F}^{2} ,

t > ik − i(i + 1)/2 remains invariant till the end of the current sweep, and that A_{k−1,k}^{[K]} is the null matrix. Hence, we obtain

S^{2} (A^{[K]}) = ∑_{i=1}^{k−2} ∑_{j=i+1}^{k} ∥A_{ij}^{[i+1,i+2]}∥_{F}^{2} ≤ S^{2} (A) ∑_{i=1}^{k−2} ∑_{j=i+2}^{k} Φ_{ij}^{2}

which proves the lemma. □

5. Cubic Convergence

In this section we prove the cubic convergence of the method P_{Θ}. First we derive basic tools for analyzing the effect of the Paardekooper steps which follow the sequences Ω_{5}, Ω_{6} and Ω_{7}.

5.1. An Auxiliary Lemma. Let G = (G_{ij}) , 1 ≤ i, j ≤ l, be the Paardekooper partition of a skew-symmetric matrix G of order 2l. We write

G = [ \begin{matrix}G_{1} & G_{3} \\ −G_{3}^{T} & G_{2} \end{matrix} ] ,

where G_{1} and G_{2} are square matrices of order 2l_{1} and 2l_{2}, respectively. Then the set

Q = {(p, q) : 1 ≤ p ≤ l_{1} , l_{1} + 1 ≤ q ≤ l_{1} + l_{2} = l} ,

is associated with the submatrix G_{3}. Note that L = l_{1}l_{2}, is the number of 2-by-2 blocks G_{ij} in G_{3} and the number of pairs in Q.

Lemma 4. If L Paardekooper steps are applied to G according to the sequence Ω_{R} (Q), then

(i) S (G_{k}^{[L]}) ≤ S (G_{k}) + \sqrt{2}S (G) \sqrt{∑_{i=1}^{l_{1}} ∑_{j=1}^{l_{2}} Φ_{i,l_{1}+j}^{2}}, k = 1, 2.

(ii) ∥G_{3}^{[L]}∥_{F} ≤ [S (G_{1}) + S (G_{2}) + 3.15S (G) \sqrt{∑_{i=1}^{l_{1}} ∑_{j=1}^{l_{2}} Φ_{i,l_{1}+j}^{2}}] \sqrt{∑_{i=1}^{l_{1}} ∑_{j=1}^{l_{2}} Φ_{i,l_{1}+j}^{2}}.

Proof. . (i) Let k = 1. By an analysis similar to that used by Wilkinson [14] together with Lemma 2, one obtains for 1 ≤ i < j ≤ l_{1}

∥G_{ij}^{[L]}∥_{F} ≤ ∥G_{ij}∥_{F} + ∑_{t=1}^{l_{2}} ∥G_{j,l_{1}+t}^{[l_{2}(i−1)]}∥_{F} Φ_{i,l_{1}+t} + ∑_{t=1}^{l_{2}} ∥G_{i,l_{1}+t}^{[l_{2}(j−1)]}∥_{F} Φ_{j,l_{1}+t}.

The argument now proceeds just as in [12]. The proof of (i) in the case k = 2 and of (ii) follows by imitating the proofs of the lemmas 3.5 and 3.6 in [12]. □

5.2. Clusters and Asymptotic Assumptions. We denote the eigenvalues of A by ±ıν_{j}, where ν_{j} ≥ 0 for 1 ≤ j ≤ k and ı = \sqrt{−1}. We shall assume that all the eigenvalues of A are simple or double. In the case of simple eigenvalues we allow for clusters of two close eigenvalues. The clusters are defined by help of the quantity δ which is computed via the following algorithm (cf. [6])

Algorithm Delta  Input date are nonnegative 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.

^{3δ} ν_{7}  ν_{6} ν_{5}  ν_{4} = ν_{3}  ν_{2}  ν_{1} 0 Figure 1: k = 7

Note that δ > 0 implies at most double eigenvalues, with the exception of the zero eigenvalue which can be at most quadruple. If δ > 0, the sets (cf. [6])

C_{i} = {ν_{j} : |ν_{j} − ν_{i}| < 3δ, 1 ≤ j ≤ k} , 1 ≤ i ≤ k,

are nonempty (since ν_{i} ∈ C_{i}), every two of them are either identical or disjoint and each of them can contain at most two (identical or distinct) eigenvalues. We call C_{i}− cluster of ν_{i} even though ν_{i} can be the only element of C_{i}. If A has simple eigenvalues λ_{i} (A) , 1 ≤ i ≤ n, then

3δ ≥ min_{i≠j} |ν_{i} − ν_{j}| ≥ min_{i≠j} |λ_{i} (A) − λ_{j} (A) |

and the orders of magnitude of the three quantities above can vary significantly (see Figure 1).

Throughout this section we let A^{[0]} = A and A^{[r]}, r ≥ 1, denotes the matrix obtained after the rth step of the method P_{Θ} defined in Section 3. For δ > 0 are well defined quantities (cf. [12])

Γ_{1}^{[r]} = ^{4} \sqrt{\frac{S^{2} (A_{1}^{[r]}) + S^{2} (A_{2}^{[r]})}{δ^{2}},} Γ_{2}^{[r]} = ^{3} \sqrt{\frac{∥A_{3}^{[r]}∥_{F}}{δ}} , r ≥ 0, (7)

and Γ^{[r]} = max {Γ_{1}^{[r]}, Γ_{2}^{[r]}} , r ≥ 0, (8)

where A_{1}^{[r]}, A_{2}^{[r]} and A_{3}^{[r]} are from the partition (5) of A^{[r]}. If r = 0 we simply write Γ_{1}, Γ_{2} and Γ. We make three asymptotic assumptions. The first one refers to the magnitudes of ∥A_{ij}∥_{F}, i < j:

A1: Γ < \frac{1}{2} .

It implicitly presumes δ > 0. The assumption A1 and the relations (7) and (8) imply

S^{2} (A) = S^{2} (A_{1}) + S^{2} (A_{2}) + ∥A_{3}∥_{F} < \frac{δ^{2}}{16} + \frac{δ^{2}}{64} = \frac{5}{64}δ^{2}. (9)

This is a stronger assumption than S (A) < δ/2 which has been used in [6]. Hence all we have proved in [6] for the Paardekooper method remains valid. In fact we can sharpen the angle bounds using the stronger condition (9). If the eigenvalue pairs ±ıν_{j} are appropriately ordered, an application of the Hoffman-Wielandt theorem (see [8, 5]) implies ||a_{2j−1,2j}| − ν_{j}| < (\sqrt{5}/8) δ, 1 ≤ j ≤ k. If ν_{j} ∈ C_{i} (ν_{j} ∉ C_{i}) we say that the block A_{jj} is (is not) associated with the block A_{ii} and write A_{jj} ∼ A_{ii} (A_{jj} ≁ A_{ii}). Since each C_{i} has at most two elements, at most two diagonal blocks can be associated. From [6, Lemma 2.1] follows the preservation of the association of diagonal blocks under Paardekooper steps.

Following the proof of [6, Lemma 2.2] and using the relation (9) one can prove^{1} that Φ_{pq} ≤ 0.85 ∥A_{pq}∥ /δ, provided that A_{pp} is not associated with A_{qq}. Hence, by the same argument as in [6, relation (3.14)] we have

\sqrt{∑_{r∈N_{s}^{t}} Φ_{p(r),q(r)}^{2}} ≤ 0.85\frac{S (A^{[s]} (t))}{δ}, 0 ≤ s < t. (10)

Here N_{s}^{t} stands for the set {r : s < r ≤ t and A_{p(r),p(r)} ≁ A_{q(r),q(r)}} and A^{[s]}(t) is the smallest principal submatrix of A^{[s]} containing all A_{p(r)q(r)}^{[s]}, s < r ≤ t. Recall that (p (r) , q (r)) is the pivot pair at step r.

The remaining two asymptotic assumptions are:

1An easy calculation shows that the condition (9) reduces the constants b, c and a in [6, Lemma 2.2]. In particular, b = 0.756 > max_{0≤x≤1}(x + 1.0561)/(2.7204915 − 0.030328\sqrt{1 − x^{2}}), c = 0.3882 and a = 0.85 > \sqrt{b^{2} + c^{2}}.

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

A3: The positive integers k_{1}, k_{2}, k_{3}, and k_{4} from the block partitions (5) are so chosen that every pair of associated diagonal 2-by-2 blocks is completely contained within one of the blocks B_{11}, B_{22}, B_{33}, or B_{44}.

To obtain sharp cubic convergence bounds we make use of the quantities Φ_{i}^{[r]},

Φ_{i}^{[r]} = \sqrt{∑_{\begin{matrix}(p,q)∈Ω_{i} \\ q≠p+1 \end{matrix}} Φ_{pq}^{2}}, 1 ≤ i ≤ 4, Φ_{i}^{[r]} = \sqrt{∑_{(p,q)∈Ω_{i}} Φ_{pq}^{2}}, 5 ≤ i ≤ 7. (11)

Here r indicates that the first of Jacobi annihilators J_{pq}, (p, q) ∈ Ω_{i} is applied to A^{[r]}. Since the association of the diagonal 2-by-2 blocks is preserved, the assumption A2 and the relation (10) imply

Φ_{i}^{[r]} ≤ 0.85\frac{S (B_{ii}^{[r]})}{δ}, 1 ≤ i ≤ 4, r ≥ 0.

Similarly, the assumption A3 and the relation (10) imply

Φ_{5}^{[r]} ≤ 0.85\frac{S (A_{1}^{[r]})}{δ}, Φ_{6}^{[r]} ≤ 0.85\frac{S (A_{2}^{[r]})}{δ}, Φ_{7}^{[r]} ≤ 0.85\frac{S (A^{[r]})}{δ}, r ≥ 0. (12)

5.3. Cubic Convergence Theorem. The central result of the paper is

Theorem 5. Let A be a skew-symmetric matrix of order n = 2k satisfying the assumptions A1, A2, and A3. Let A^{[0]} = A, A^{[1]}, …, A^{[M]} be obtained by applying the method P_{Θ} to A. Then

Γ_{1}^{[M]} < 1.963 Γ^{3} and Γ_{2}^{[M]} < 1.775 Γ^{3}.

Proof. . Recall from Definition 1 that one quasi-cycle of the method P_{Θ} is defined by the sequence Θ. It is convenient to break Θ into four parts

Θ = Θ_{1}, Θ_{2}, Θ_{3}, Θ_{4}

where

Θ_{1} = Ω_{1}, Ω_{2}, Ω_{3}, Ω_{4}, Θ_{2} = Ω_{5}, Ω_{6}, Θ_{3} = Θ_{1}, Θ_{4} = Ω_{7}.

With each sequence Θ_{i}, 1 ≤ i ≤ 4, we associate a stage of the process. By μ_{i} we denote the cardinality of the set P_{i}, 1 ≤ i ≤ 7, from Section 3.

Stage 1: Here we consider the part of the process defined by the sequence Θ_{1} = Ω_{1}, Ω_{2}, Ω_{3}, Ω_{4}. We first concentrate on the sequence Ω_{1}. Applying Lemma 3 to B_{11} and using the relation (11) with i = 1 we obtain

S (B_{11}^{[μ_{1}]}) ≤ S (B_{11}) Φ_{1}^{[0]}.

Since B_{11}^{[μ_{1}]} is not affected by Jacobi annihilators which follow the sequence Ω_{2}, Ω_{3}, Ω_{4}, we can replace B_{11}^{[μ_{1}]} by B_{11}^{[r_{1}]}, where r_{1} = μ_{1} + μ_{2} + μ_{3} + μ_{4}. In a similar way we obtain

S (B_{ii}^{[r_{1}]}) ≤ S (B_{ii}) Φ_{i}^{[∑_{j=1}^{i−1} μ_{j}]}, 1 ≤ i ≤ 4, (13) where empty sum is considered zero. Note that

S (A_{k}^{[r_{1}]}) ≤ S (A_{k}) k = 1, 2. (14)

Stage 2: In this stage we consider the part of the process defined by the sequence Θ_{2} = Ω_{5}, Ω_{6}. We first concentrate on the sequence Ω_{5}. Applying Lemma 4 (i) to G = A_{1}^{[r_{1}]} with G_{1} = B_{11}^{[r_{1}]}, G_{2} = B_{22}^{[r_{1}]} and G_{3} = B_{12}^{[r_{1}]}, we obtain

S (B_{ii}^{[r_{1}+μ_{5}]}) ≤ S (B_{ii}^{[r_{1}]}) + \sqrt{2}S (A_{1}^{[r_{1}]}) Φ_{5}^{[r_{1}]}, i = 1, 2.

Because B_{ii}^{[r_{1}+μ_{5}]}, i = 1, 2 are not affected by Jacobi annihilators which follow the sequence Ω_{6}, we can replace the superscript r_{1} + μ_{5} by r_{2}, where r_{2} = r_{1} + μ_{5} + μ_{6}. Using the relations (13) and (14) we obtain

S (B_{ii}^{[r_{2}]}) ≤ S (B_{ii}) Φ_{i}^{[∑_{j=1}^{i−1} μ_{j}]} + \sqrt{2} S (A_{1}) Φ_{5}^{[r_{1}]}, i = 1, 2.

Note that S (B_{ii}) ≤ S (A_{1}) , i = 1, 2. Using the Cauchy-Schwarz inequality and the relation (10) we obtain

(15)  max \{S (B_{11}^{[r_{2}]}) , S (B_{22}^{[r_{2}]})\} ≤ \sqrt{3} S (A_{1}) \sqrt{(max \{Φ_{1}^{[0]}, Φ_{2}^{[μ_{1}]}\})^{2} + (Φ_{5}^{[r_{1}]})^{2}}
≤ \sqrt{3} S (A_{1}) · 0.85 \frac{S (A_{1})}{δ} ≤ 1.4723 \frac{S^{2} (A_{1})}{δ}.

Similarly, we obtain

max \{S (B_{33}^{[r_{2}]}) , S (B_{44}^{[r_{2}]})\} ≤ 1.4723 \frac{S^{2} (A_{2})}{δ}.

Using Lemma 4 (ii) with G = A_{1}^{[r_{1}]}, G_{1} = B_{11}^{[r_{1}]}, G_{2} = B_{22}^{[r_{1}]} and G_{3} = B_{12}^{[r_{1}]}, we obtain

∥B_{12}^{[r_{2}]}∥_{F} ≤ [S (B_{11}^{[r_{1}]}) + S (B_{22}^{[r_{1}]}) + 3.15 S (A_{1}^{[r_{1}]}) Φ_{5}^{[r_{1}]}] Φ_{5}^{[r_{1}]}.

Here we have used the fact that B_{12}^{[r_{2}]} = B_{12}^{[r_{1}+μ_{5}]}. Using the relations (13) , (14), the Cauchy-Schwarz inequality, and the relation (10) , we obtain

∥B_{12}^{[r_{2}]}∥_{F} ≤ [S (B_{11}) Φ_{1}^{[0]} + S (B_{22}) Φ_{2}^{[μ_{1}]} + 3.15 S (A_{1}) Φ_{5}^{[r_{1}]}] Φ_{5}^{[r_{1}]} (16)

≤ \sqrt{S^{2} (B_{11}) + S^{2} (B_{22}) + 3.15^{2} S^{2} (A_{1})}\sqrt{(Φ_{1}^{[0]})^{2} + (Φ_{2}^{[μ_{1}]})^{2} + (Φ_{5}^{[r_{1}]})^{2}} Φ_{5}^{[r_{1}]}
≤ \sqrt{1 + 3.15^{2}}S (A_{1}) · 0.85 · \frac{S (A_{1})}{δ} Φ_{5}^{[r_{1}]} ≤ 2.81 \frac{S^{2} (A_{1})}{δ} Φ_{5}^{[r_{1}]}.

Similarly, we have

∥B_{34}^{[r_{2}]}∥_{F} ≤ 2.81 \frac{S^{2} (A_{2})}{δ} Φ_{6}^{[r_{1}+μ_{5}]}.

Stage 3: Since Θ_{3} = Θ_{1}, this stage is entirely analogous to stage 1. Let r_{3} = r_{1} + r_{2}. Using Lemma 3 and the relation (11) we obtain

S (B_{ii}^{[r_{3}]}) ≤ S (B_{ii}^{[r_{2}]}) Φ_{i}^{[r_{2}+∑_{j=1}^{i−1} μ_{j}]}, 1 ≤ i ≤ 4. (17)

Note that ∥B_{12}^{[r_{3}]}∥_{F} = ∥B_{12}^{[r_{2}]}∥_{F}. Using the relations (17) , (16), (15) , (10) and (14) , we obtain

(18)  S^{2} (A_{1}^{[r_{3}]}) = S^{2} (B_{11}^{[r_{3}]}) + S^{2} (B_{22}^{[r_{3}]}) + ∥B_{12}^{[r_{2}]}∥_{F}^{2}
≤ S^{2} (B_{11}^{[r_{2}]}) (Φ_{1}^{[r_{2}]})^{2} + S^{2} (B_{22}^{[r_{2}]}) (Φ_{2}^{[r_{2}+μ_{1}]})^{2} + 2.81^{2} \frac{S^{4} (A_{1})}{δ^{2}} (Φ_{5}^{[r_{1}]})^{2}
≤ 1.4723^{2} \frac{S^{4} (A_{1})}{δ^{2}} [(Φ_{1}^{[r_{2}]})^{2} + (Φ_{2}^{[r_{2}+μ_{1}]})^{2}] + 2.81^{2} \frac{S^{4} (A_{1})}{δ^{2}} (Φ_{5}^{[r_{1}]})^{2}
≤ 2.81^{2} \frac{S^{4} (A_{1})}{δ^{2}} [(Φ_{1}^{[r_{2}]})^{2} + (Φ_{2}^{[r_{2}+μ_{1}]})^{2} + (Φ_{5}^{[r_{1}]})^{2}]

≤ (2.81 · 0.85)^{2} \frac{S^{6} (A_{1}^{[r_{1}]})}{δ^{4}} ≤ 5.705 \frac{S^{6} (A_{1})}{δ^{4}}.

In a similar way, we obtain

S^{2} (A_{2}^{[r_{3}]}) ≤ 5.705 \frac{S^{6} (A_{2})}{δ^{4}}. (19)

The relations (18) , (19) and (7) yield

S^{2} (A_{1}^{[r_{3}]}) + S^{2} (A_{2}^{[r_{3}]}) ≤ 5.705 [\frac{S^{2} (A_{1}) + S^{2} (A_{2})}{δ^{2}}]^{3} δ^{2} ≤ 5.705 Γ_{1}^{12}δ^{2}. (20)

Since ∥A_{3}^{[r_{3}]}∥_{F} = ∥A_{3}∥_{F}, the relations (20) , (7) , (8) and the assumption A1 imply

(21)  S^{2} (A^{[r_{3}]}) = S^{2} (A_{1}^{[r_{3}]}) + S^{2} (A_{2}^{[r_{3}]}) + ∥A_{3}∥_{F}^{2} ≤ 5.705 Γ_{1}^{12}δ^{2} + Γ_{2}^{6}δ^{2}
≤ \frac{5.705}{2^{6}}Γ_{1}^{6}δ^{2} + Γ_{2}^{6}δ^{2} ≤ 1.09 Γ^{6}δ^{2}.

Stage 4: In this stage, which is similar to stage 2, we consider the remaining μ_{7} steps defined by the sequence Θ_{4} = Ω_{7}. Using Lemma 4 (i) with G = A^{[r_{3}]}, G_{1} = A_{1}^{[r_{3}]}, G_{2} = A_{2}^{[r_{3}]}, G_{3} = A_{3}^{[r_{3}]} and the relations (11) and (12), we obtain

(22)  S (A_{k}^{[M]}) ≤ S (A_{k}^{[r_{3}]}) + \sqrt{2} S (A^{[r_{3}]}) Φ_{7}^{[r_{3}]}
≤ S (A_{k}^{[r_{3}]}) + \sqrt{2} · 0.85 \frac{S^{2} (A^{[r_{3}]})}{δ}, k = 1, 2.

Using the relations (7), (22), the Cauchy-Schwarz inequality and the relations (20), (21) and (8) we obtain

(Γ_{1}^{[M]})^{4} = \frac{S^{2} (A_{1}^{[M]}) + S^{2} (A_{2}^{[M]})}{δ^{2}} ≤ 2\frac{S^{2} (A_{1}^{[r_{3}]}) + S^{2} (A_{2}^{[r_{3}]})}{δ^{2}} + 1.7^{2} \frac{S^{4} (A^{[r_{3}]})}{δ^{4}}
≤ 2 · 5.705Γ_{1}^{12} + 1.7^{2} · 1.09^{2}Γ^{12} ≤ 14.844Γ^{12}.

The last inequality implies Γ_{1}^{[M]} < 1.963 Γ^{3}.

To prove the second assertion of the theorem we apply Lemma 4 (ii) to G = A^{[r_{3}]} with G_{1} = A_{1}^{[r_{3}]}, G_{2} = A_{2}^{[r_{3}]} and G_{3} = A_{3}^{[r_{3}]}. Using the relation (11) with i = 7 , the Cauchy-Schwarz inequality, and the relations (12) , (20) , and (21) , we obtain

∥A_{3}^{[M]}∥_{F} ≤ [S (A_{1}^{[r_{3}]}) + S (A_{1}^{[r_{3}]}) + 3.15 S (A^{[r_{3}]}) Φ_{7}^{[r_{3}]}] Φ_{7}^{[r_{3}]} ≤
≤ [\sqrt{2} \sqrt{S^{2} (A_{1}^{[r_{3}]}) + S^{2} (A_{1}^{[r_{3}]})} + 3.15 S (A^{[r_{3}]}) · 0.85\frac{S (A^{[r_{3}]})}{δ}] · 0.85\frac{S (A^{[r_{3}]})}{δ}
≤ [\sqrt{2} \sqrt{5.705} Γ_{1}^{6}δ + 3.15 · 0.85 · 1.09 Γ^{6}δ] 0.85\sqrt{1.09}Γ^{3} ≤ 5.588 Γ^{9}δ.

Using the relation (7) , we obtain the second assertion of the theorem

Γ_{2}^{[M]} ≤ ^{3} \sqrt{5.588} Γ^{3} < 1.775 Γ^{3}. □

From the relations (7) and (8) and the assumption A1 we see that

Γ^{3} ≤ \sqrt{Γ_{1}^{4} + Γ_{2}^{6}} = \frac{S (A)}{δ} ≤ \frac{\sqrt{5}}{2}Γ^{2}.

By Theorem 5, Γ reduces cubically per quasi-sweep. Hence S (A) /δ has upper and lower bounds which converge cubically per quasi-sweep. This implies the ultimate cubic convergence of the method P_{Θ}.

Corollary 6. Let the assumptions of Theorem 5 hold.

(i) If Γ = Γ_{1}, then \frac{S (A^{[M]})}{δ} ≤ 3.92 [\frac{S (A)}{δ}]^{3}

(ii) If Γ = Γ_{2}, then \frac{S (A^{[M]})}{δ} ≤ 3.92 [\frac{S (A)}{δ}]^{2}

(iii) If Γ_{2} = γΓ_{1}, then \frac{S (A^{[M]})}{δ} ≤ 3.92 max \{1, γ^{6}\} [\frac{S (A)}{δ}]^{3} .

Proof. . By the relations (7) , (8) , Theorem 5 and the assumption A1, we have

(23)  \frac{S^{2} (A^{[M]})}{δ^{2}} = (Γ_{1}^{[M]})^{4} + (Γ_{2}^{[M]})^{6} ≤ (1.963 Γ^{3})^{4} + (1.775Γ^{3})^{6}
≤ 1.963^{4} Γ^{12} + 1.775^{6} · \frac{1}{2^{6}} Γ^{12} ≤ 15.338 Γ^{12}.

(i) If Γ = Γ_{1}, then

Γ^{6} = (Γ^{4})^{\frac{3}{2}} = (Γ_{1}^{4})^{\frac{3}{2}} = [\frac{S^{2} (A_{1}) + S^{2} (A_{2})}{δ^{2}}]^{\frac{3}{2}} ≤ [\frac{S^{2} (A)}{δ^{2}}]^{\frac{3}{2}} = \frac{S^{3} (A)}{δ^{3}},

and the result follows by the relation (23) . (ii) If Γ = Γ_{2}, then Γ^{6} = Γ_{2}^{6} = \frac{∥A_{3}∥_{F}^{2}}{δ^{2}} ≤ \frac{S^{2}(A)}{δ^{2}}, and the result follows again by the relation (23) . (iii) Since Γ = max {1, γ} Γ_{1}, the proof is quite similar to that of (i). □

From Corollary 6, if Γ = Γ_{2} and γ is not of moderate size the cubic reduction of S (A) /δ per quasi-sweep may fail to occur. However, cubic convergence of Γ actually forces cubic convergence of S (A) /δ. Once the asymptotic assumptions are fulfilled, the quadratic reduction of S (A) /δ can occur just one time. The argument is same as that presented in [12].

It is easy to see that the number of Jacobi annihilators in one quasi-sweep of the method P_{Θ} is minimized when k_{1}, k_{2}, k_{3}, and k_{4} defined in section 3 are approximately the same number. In such a case, one quasi-sweep is approximately equal to 1.25 of ordinary sweeps.

5.4. Numerical Examples. We have made numerical tests for the method P_{Θ}. All the computations were done on VAX6000 Model 500 using quadruple precision arithmetic with machine precision eps ≃ 9.63 x 10^{−35}. The results have confirmed all the conclusions we stated above. Here are some typical examples.

Example 7. Let n = 32, A = (a_{ij}) , A = −A^{T} , a_{2i−1,2i} = i, 1 ≤ i ≤ 16 and a_{ij} = 10^{−8} for all other elements in the upper triangular part of A. From the computed eigenvalues we have found out that δ = 0.333333333333333339318783068783067 to the arithmetic accuracy. We chose k_{1} = k_{2} = k_{3} = k_{4} = 4. Applying P_{Θ} to A we have computed Γ^{[tM]}, S (A^{[tM]}) /δ, γ_{t} = Γ_{2}^{[tM]}/Γ_{1}^{[tM]}, and

g_{t} = \frac{Γ^{[tM]}}{1.963 (Γ^{[(t−1)M]})^{3}} , s_{t} = \frac{S (A^{[tM]}) /δ}{3.92 max \{1, (γ_{t−1})^{6}\} [S (A^{[(t−1)M]}) /δ]^{3}} ,

after every quasi-sweep (t counts quasi-sweeps). From the first row of the table and from the fact that a_{2i−1,2i} = i, 1 ≤ i ≤ 16, we see that the assumptions A1, A2, and A3 are satisfied. Hence by Theorem 5 and Corollary 6 (iii), g_{t} and s_{t} cannot be greater than 1. Applying P_{Θ} to A we have obtained the following results:

t  Γ^{[tM]}  g_{t}  γ_{t}  S (A^{[tM]}) /δ  s_{t}
0 7.830 · 10^{−3}  11.685 6.573 · 10^{−7}
1 8.736 · 10^{−8} 9.272 · 10^{−2} 0.4350 7.632 · 10^{−15} 2.694 · 10^{−3}
2 1.727 · 10^{−23} 1.323 · 10^{−2} 0.4250 2.982 · 10^{−46} 1.711 · 10^{−4}
3 1.627 · 10^{−70} 1.610 · 10^{−2} 0.3142 2.649 · 10^{−140} 2.550 · 10^{−4}

We see that in the case Γ = Γ_{2} we cannot guarantee cubic reduction of S (A) /δ per quasi-sweep. But as we mentioned eariler, the quadratic reduction of S (A) /δ can occur just one time. □

In the following example we consider double eigenvalues.

Example 8. We have generated a matrix of order n = 33 from the real Schur form whose nontrivial entries are 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8 by a series of about 2.5 n^{2} similarity transformations with plane rotations. For each matrix generation we have used the row-cyclic ordering and the rotation angles were chosen randomly from the interval [−1/5000, 1/5000] . Such angle choices ensured that A satisfies the assumptions A1 and A2. Obviously δ = 1/3. To ensure that A satisfies the assumption A3, we have chosen k_{1} = k_{2} = k_{3} = 4, and k_{4} = 5. Applying P_{Θ} to A we have obtained the following results displayed in the table below:

t  Γ^{[tM]}  g_{t}  γ_{t}  S (A^{[tM]}) /δ  s_{t}
0 4.458 · 10^{−1}  1.5374 1.221 · 10^{−1}
1 9.226 · 10^{−3} 5.305 · 10^{−2} 0.4734 8.511 · 10^{−5} 9.025 · 10^{−4}
2 1.789 · 10^{−8} 1.160 · 10^{−2} 0.4809 3.200 · 10^{−16} 1.324 · 10^{−4}
3 3.046 · 10^{−25} 2.710 · 10^{−2} 0.2804 9.278 · 10^{−50} 7.222 · 10^{−4}
4 1.767 · 10^{−75} 3.186 · 10^{−2} 0.1475 3.124 · 10^{−150} 9.973 · 10^{−4}

Note the cubic reduction of Γ and of S (A) /δ. A similar behavior can be seen in the presence of clusters of two simple eigenvalues. □

Next we show that the assumption A3 cannot be dropped.

Example 9. We have used the matrix from Example 8. But we have chosen k_{1} = 3, k_{2} = 6, k_{3} = k_{4} = 4 so that the assumption A3 are violated even though the assumptions A1 and A2 are still true. Applying P_{Θ} to A we have obtained the following results:

t  Γ^{[tM]}  g_{t}  γ_{t}  S (A^{[tM]}) /δ  s_{t}
0 4.508 · 10^{−1}  1.586 · 10^{0} 1.221 · 10^{−1}
1 9.429 · 10^{−3} 5.243 · 10^{−2} 4.752 · 10^{−1} 8.891 · 10^{−5} 7.808 · 10^{−4}
2 1.886 · 10^{−8} 1.146 · 10^{−2} 4.975 · 10^{−1} 3.558 · 10^{−16} 1.292 · 10^{−4}
3 8.873 · 10^{−15} 6.735 · 10^{8} 1.951 · 10^{−4} 7.872 · 10^{−29} 4.460 · 10^{17}
4 4.301 · 10^{−35} 3.138 · 10^{7} 3.461 · 10^{−8} 1.851 · 10^{−69} 9.685 · 10^{14}
5 3.125 · 10^{−107} 1.998 · 10^{−4} 2.976 · 10^{−9} 9.766 · 10^{−214} 3.925 · 10^{−8}

Note that Γ^{[tM]} is by no means cubic reduction of Γ^{[(t−1)M]} for t = 3, 4. Also note that S (A^{[tM]}) /δ is not cubic reduction of S (A^{[(t−1)M]}) /δ even though Γ^{[(t−1)M]} = Γ_{1}^{[(t−1)M]}, t = 3, 4. Hence Theorem 5 and Corollary 6 cannot be true without the assumption A3. □

Finally, we show that cubic convergence can fail to hold when multiplicities of the eigenvalues are greater that two. This is a sharp contrast to the case of symmetric Jacobi method under the same quasi-cyclic strategy. If 3 δ̃ is the minimum distance between two distinct eigenvalues of a symmetric matrix A, then under the appropriate assumptions the quasi-cyclic Jacobi method yields cubic reduction of S (A) /δ̃ [12].

Example 10. By the same procedure as in Example 8 we have generated A, of order n = 24, from the Murnaghan form whose nontrivial entries are 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4. If we set 3 δ̃ = min_{ν_{i}≠ν_{j}} |ν_{i}−ν_{j}|, then we have δ̃ = 1/3. We chose k_{1} = k_{2} = k_{3} = k_{4} = 6. Let Γ̃^{[tM]}, g̃_{t}, γ̃_{t}, and s̃_{t} be defined using δ̃ in place of δ. It was checked that the initial matrix A satisfies the assumptions A1, A2, and A3 using Γ̃ and δ̃ in place of Γ and δ, respectively. Applying P_{Θ} to A we have obtained the following results displayed in the table below:

t  Γ̃^{[tM]}  g̃_{t}  γ̃_{t}  S (A^{[tM]}) /δ̃  s̃_{t}
0 3.214 · 10^{−1}  1.916 · 10^{0} 4.351 · 10^{−2}
1 7.021 · 10^{−3} 1.077 · 10^{−1} 4.828 · 10^{−1} 4.929 · 10^{−5} 3.085 · 10^{−3}
2 4.301 · 10^{−5} 6.331 · 10^{1} 3.977 · 10^{−2} 1.850 · 10^{−9} 3.940 · 10^{3}
3 9.336 · 10^{−9} 5.979 · 10^{4} 2.334 · 10^{−4} 8.716 · 10^{−17} 3.514 · 10^{9}
4 1.503 · 10^{−12} 9.412 · 10^{11} 5.023 · 10^{−9} 2.260 · 10^{−24} 8.708 · 10^{23}

We note that S (A) /δ̃ fails to converge cubically (quadratically). A similar behavior can be seen in the presence of clusters with more than two simple eigenvalues. This behavior is caused by a special structure of a skew-symmetric matrix which is almost in Murnaghan form. This structure is derived in [3], but a more elaborated discussion can be found in [7]. □

The last example shows that δ defined by Algorithm Delta is the optimal (that is, the largest) one under which the cubic convergence of the method P_{Θ} is guaranteed.

References

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

[2] D. Hacon, Jacobi’s Method for Skew-Symmetric Matrices, SIAM J. on Matrix Analysis and Applications 14, 619–628, 1993.

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

[4] V. Hari, On Sharp Quadratic Convergence Bounds for the Serial Jacobi Methods, Numer. Math. 60, 375-406, 1991.

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

[6] V. Hari and N. H. Rhee, On the Quadratic Convergence of the Paardekooper Method II, Glasnik Matematički 27 (47), 369–391, 1992.

[7] V. Hari and N. H. Rhee, A Matrix Pair of an Almost Dioagonal Skew-Symmetric Matrix and a Symmetric Positive Definite Matrix, Linear Algebra and Its Applications 195, 117–144, 1993.

[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] W. F. Mascarenhas, On the Convergence of the Jacobi Method for Arbitrary Orderings, PhD thesis, M.I.T., 1991.

[10] F. D. Murnaghan and A. Winter, A Canonical Form for Real Matrices under Orthogonal Transformations, USA Nat. Acad. Sci. 17, 417–420, 1931.

[11] M. H. C. Paardekooper, An Eigenvalue Algorithm for Skew-Symmetric Matrices, Numer. Math. 17, 189–202, 1971.

[12] N. H. Rhee and V. Hari, On the Global and Cubic Convergence of a Quasi-Cyclic Jacobi Method, Numer. Math. 66, 97–122, 1993.

[13] 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.

[14] J. H. Wilkinson, Note on the Quadratic Convergence of the Cyclic Jacobi Process, Numer. Math. 4, 296–300, 1962.