Quadratic Convergence of Scaled Iterates by Kogbetliantz Method

Vjeran Hari^{∗} and Josip Matejaš^{†}

Abstract. A sharp quadratic convergence bound of scaled iterates for the serial Kogbetliantz method is derived. Iterates are symmetrically scaled by diagonal matrices so that diagonal elements are ones. The result is obtained for a scaled diagonally dominant complex triangular matrix with multiple singular values. The estimate depends on the relative separation of the singular values.

Key words. Kogbetliantz method, scaled matrices, quadratic convergence.

AMS subject classification. 65F15, 65G05.

0 Introduction

Kogbetliantz method for computing SVD of triangular matrices is a reliable simple method well suited for ordinary [9, 3] and parallel [9, 13] computers. Our numerical tests give directions to conclusion that it is relatively accurate at least for matrices of the form DA, AD and DAD where D is diagonal and A is well conditioned. Using recent estimates [19] on the relative distance between diagonal elements of a scaled diagonally dominant triangular matrix and the corresponding singular values, one can conclude that the off-norm of the symmetrically scaled iterate is a good measure for advancing of the algorithm. The phenomenon is quite similar to the known one (see [4]) for the standard Jacobi method for Hermitian positive definite matrices. The stopping criterion based on the new measure can ensure the relative accuracy of the computed singular values and sometimes (e.g. graded matrices) it can terminate the process sooner than the classical stopping criterion.

These observations have motivated a research of scaled iterates by diagonalization methods (see [16, 17, 18, 20]) and this paper is part of that program. In [20] we have obtained usable estimates of scaled iterates by Kogbetliantz method and have made a sketch of the quadratic convergence proof in the case of simple singular values. Here we provide a detailed quadratic convergence proof in the general situation of multiple singular values. The proof is valid for the case of a complex triangular initial matrix and the complex Kogbetliantz algorithm. Here we use the technique from [18, 16], hence the reader is suggested to be acquainted with [17, 18] (or [16]) and [20].

∗Department of Mathematics, University of Zagreb, Bijenička 30, 10000 Zagreb, CROATIA

†Faculty of Economics, University of Zagreb, Kennedyjev trg 6, 10000 Zagreb, CROATIA

The paper is organized as follows. In Section 1 we present recent results on the structure of scaled diagonally dominant triangular matrices which use relative gaps in the set of singular values. In Section 2 we shortly describe the Kogbetliantz method and in Section 3 we provide some preparatory results. In Section 4 we formulate the asymptotic assumptions and provide more auxiliary results. In Section 5 we prove the main result, Theorem 5.1 and in Section 6 we provide some numerical examples.

1 Structure of scaled almost diagonal matrices

Let C^{n×n} denote the set of complex n by n matrices and let B ∈C^{n×n}. Let

σ_{1} = σ_{2} = ⋯ = σ_{s_{1}} > σ_{s_{1}+1} = ⋯ = σ_{s_{2}} > ⋯ > σ_{s_{p−1}+1} = … = σ_{s_{p}} (1.1)

be the singular values of B. Then p is the number of distinct singular values of B and for each i, 1 ≤ i ≤ p, n_{i} = s_{i} − s_{i−1} (s_{0} = 0, s_{p} = n) is the multiplicity of σ_{s_{i}}. Let us define the sets of indices

N_{i} = {t ∈ N : s_{i−1} + 1 ≤ t ≤ s_{i}}, 1 ≤ i ≤ p .  (1.2)

According to the partition n = n_{1} + ⋯ + n_{p} we define for every X ∈C^{n×n} the block partition

X = [ \begin{matrix}X_{11} & … & X_{1p} \\ ⋮ & & ⋮ \\ X_{p1} & … & X_{pp} \end{matrix} ] , X_{ij} ∈ C^{n_{i}×n_{j}}, 1 ≤ i, j ≤ p.

In the sequel we use the following notation: Ω(X) = X−diag(X) is the offdiagonal part of X, π(X) =diag(X_{11}, …, X_{pp}) is the block-diagonal and τ(X) = X − π(X) is the off-block-diagonal part of X. The i’th diagonal block and the the corresponding off-diagonal block-row of X are denoted by

π_{i}(X) = X_{ii}, τ_{i}(X) = [X_{i1} … X_{i,i−1} X_{i,i+1} … X_{ip}], 1 ≤ i ≤ p .

By ||X||_{2} and ||X|| are denoted the spectral and the Frobenius norm of X respectively. For B with nonzero diagonal, the matrix

B_{S} = |diag(B)|^{−1/2}B|diag(B)|^{−1/2}  (1.3)

is referred to as the scaled matrix (we use here only symmetric scaling). Note that all diagonal elements of B_{S} have absolute value one.

If ||Ω(B_{S})|| ≤ α with α < 1 then B is α-scaled diagonally dominant (α-s.d.d.) matrix with respect to norm || · || (see [1]).

The relative gap function rg of two real arguments is defined by (see [12])

rg(a, b) = \{ \begin{matrix}\frac{|a − b|}{|a| + |b|}, & |a| + |b| > 0 \\ 0, & a = b = 0. \end{matrix}  (1.4)

Relative gaps in the set of singular values of B are numbers

θ_{i} = min_{j≠i^{1≤j≤p}} rg(σ_{s_{i}}, σ_{s_{j}}), 1 ≤ i ≤ p,  (1.5)

and the minimum relative gap is

θ = min_{1≤i≤p} θ_{i} .  (1.6)

Here we have used (1.1). Note that θ_{i} ≤ 1 for all 1 ≤ i ≤ p.

The structure of an α-s.d.d. matrix with multiple singular values is given in [19]. Here we consider only triangular and essentially triangular matrices (satisfying b_{rt}b_{tr} = 0, r ≠ t). Therefore, we recall [19, Corollary 6] which reveals the special structure of an α-s.d.d. essentially triangular matrix^{1}.

Theorem 1.1 Let B ∈C^{n×n} be essentially triangular matrix with positive diagonal elements. Let

(a)  (∀q ∈ N_{i}) (∀r ∈ N_{j}) i < j ⇒ b_{qq} > b_{rr} ,
(b)  ||Ω(B_{S})||_{2} < \frac{θ}{3},

where N_{i}, B_{S} and θ are defined by the relations (1.2), (1.3) and (1.6), respectively. Then

(i)  ∑_{t=s_{i−1}+1}^{s_{i}} |1 − \frac{σ_{s_{i}}}{b_{tt}} |^{2} + \frac{1}{2}||Ω(π_{i}(B_{S}))||_{F}^{2} ≤ \frac{4}{θ_{i}^{2}} [||Ω(π_{i}(B_{S}))||_{F}^{2}

+ ||τ_{i}(B_{S})||_{F}^{2} + ||τ_{i}(B_{S}^{∗})||_{F}^{2} ]^{2} , 1 ≤ i ≤ p,

(ii)  ∑_{t=1}^{n} |1 − \frac{σ_{t}}{b_{tt}} |^{2} + \frac{1}{2}||Ω(π(B_{S}))||_{F}^{2}

≤ \frac{4}{θ^{2}} [||Ω(B_{S})||_{F}^{2} + ||τ(B_{S})||_{F}^{2} ] · ||Ω(B_{S})||_{F}^{2} .

We see that for an α-s.d.d. matrix with α < θ/3, the off-diagonal elements of diagonal block π_{i}(B_{S}) (which is affiliated with σ_{s_{i}}) are quadratically small with respect to the norm of the off-block-diagonal part of the matrix, ||τ(B_{S})||_{F}. This structure has impact on the rate of convergence of scaled iterates by the Kogbetliantz method.

1The original assumption b_{11} ≥ ⋯ ≥ b_{nn} is replaced here by the weaker one (a).

2 Kogbetliantz method

Here we give a brief description of the serial Kogbetliantz method (see [14, 15, 7, 22] for complex triangular matrices (cf. [9, 10]) with positive diagonal elements. Kogbetliantz method for complex triangular matrices has a nice property to maintain the reality (positivity) of the initially real (positive) diagonal elements. Since we consider scaled iterates, we shall assume positive diagonal elements which satisfy the condition (a) of Theorem 1.1.

In general, nonnegative diagonal elements can easily be obtained by premultiplying a complex triangular matrix with a unitary diagonal matrix. Note that, as part of a SVD computation, QR factorization with column pivoting is often used as a pre-processing step (see also [5]). And it leads to a triangular matrix with real diagonal satisfying (up to signs) the condition (a) of Theorem 1.1. Note also that zero singular values of a triangular matrix can be deflated using Givens rotations.

The method we work with is actually the one from [20], but for the completeness of this paper we shortly describe it here.

Starting from a triangular matrix B = (b_{lm}), Kogbetliantz method generates a sequence of matrices B^{(k)} = (b_{lm}^{(k)}), k ≥ 0, by the rule

B^{(k+1)} = (U^{(k)})^{∗}B^{(k)}V^{(k)}, k ≥ 0 ; B^{(0)} = B ,  (2.1)

where U^{(k)} and V^{(k)}, k ≥ 0 are unitary plane matrices whose essential four elements are

(U^{(k)})_{ii} = (U^{(k)})_{jj} = cos φ^{(k)}, (U^{(k)})_{ij} = −\overline{(U^{(k)})_{ji}} = e^{ı ω_{k}} sin φ^{(k)},
(V^{(k)})_{ii} = (V^{(k)})_{jj} = cos ψ^{(k)}, (V^{(k)})_{ij} = −\overline{(V^{(k)})_{ji}} = e^{ı ϑ_{k}} sin ψ^{(k)}.

For all other elements, we have (U^{(k)})_{lm} = (V^{(k)})_{lm} = δ_{lm}. Here, ı is the imaginary unit, \overline{z} is the complex conjugate z and δ_{lm} is the Kronecker’s delta. The 2 × 2 matrices obtained on the intersections of rows and columns i and j of B^{(k)}, U^{(k)}, V^{(k)} are called pivot submatrices of these matrices. The pair of indices (i, j) = (i(k), j(k)) is called pivot pair and b_{ij}^{(k)} is the pivot element. The equality in (2.1) defines the k’th step or the k’th iteration of the method. In the sequel, i and j will be reserved for pivot indices.

For b_{ij}^{(k)} ≠ 0 the angles φ^{(k)}, ψ^{(k)}, ω_{k} and ϑ_{k} are determinated by the requirement

b_{ij}^{(k+1)} = b_{ji}^{(k+1)} = 0 ,  (2.2)

which ensures that ||Ω(B^{(k+1)})||^{2} = ||Ω(B^{(k)})||^{2} − |b_{ij}^{(k)}|^{2}, k ≥ 0 , holds. The relation (2.2) yields several sets of formulas for required angles, depending on whether U^{(k)} or V^{(k)} is calculated first (see [7, 22, 3, 10, 11]). Here we consider only the serial, that is, the column- and the row-cyclic pivot strategies. In the first (second) case, the off-diagonal elements are annihilated in column- (row-) wise fashion. Serial Kogbetliantz methods do not essentially ruin the triangular form. If B is upper-triangular, then within an odd (even) cycle the nontrivial off-diagonal elements push downwards (upwards) so that after each odd (even) cycle the current matrix is lower- (upper-) triangular. Similar situation is when B is lower-triangular. Within every cycle each current matrix is permutationally similar to a matrix in triangular form (PST), hence it is also essentially triangular (see [9]). Therefore, the whole process can be performed on an upper-triangular array (see [10, 11]). At step k, we can associate the uppertriangular matrix G^{(k)} = (g_{lm}^{(k)}) with the content of that array. Then, we have G^{(k)} + G^{(k)T} = B^{(k)} + B^{(k)T} , k ≥ 0 . In the following, as we mentioned earlier, we assume positive diagonal elements. We shall express the angle formulas in terms of the elements of G^{(k)}. If the right-hand transformation V^{(k)} is computed first, the relation (2.2) implies

e^{ı ϑ_{k}} tan 2ψ^{(k)} = \frac{2g_{ii}^{(k)}g_{ij}^{(k)}}{(g_{jj}^{(k)})^{2} − (g_{ii}^{(k)})^{2} + |g_{ij}^{(k)}|^{2}} ,

e^{ı ω_{k}} tan φ^{(k)} = e^{ıϑ_{k}} \frac{|g_{ij}^{(k)}| + g_{ii}^{(k)} tan ψ^{(k)}}{g_{jj}^{(k)}} = \frac{g_{jj}^{(k)}e^{ıϑ_{k}} tan ψ^{(k)}}{g_{ii}^{(k)} − |g_{ij}^{(k)}| tan ψ^{(k)}} .

Hence, we obtain the following angle formulas (cf. [22, 9, 10, 11]) for the first sweep:

\begin{matrix}e^{ı ϑ_{k}} & = & \frac{g_{ij}^{(k)}}{|g_{ij}^{(k)}|}, i.e. ϑ_{k} = arg(g_{ij}^{(k)}) \\ tan 2ψ^{(k)} & = & \frac{2g_{ii}^{(k)}|g_{ij}^{(k)}|}{(g_{jj}^{(k)})^{2}−(g_{ii}^{(k)})^{2}+|g_{ij}^{(k)}|^{2}} , ψ^{(k)} ∈ [−\frac{π}{4} , \frac{π}{4} ] , \\ z_{k} & = & \frac{g_{ij}^{(k)}+g_{ii}^{(k)}e^{ıϑ_{k}} tan ψ^{(k)}}{g_{jj}^{(k)}} = e^{ıϑ_{k}} \frac{|g_{ij}^{(k)}|+g_{ii}^{(k)} tan ψ^{(k)}}{g_{jj}^{(k)}} \\ & = & \frac{g_{jj}^{(k)}e^{ıϑ_{k}} tan ψ^{(k)}}{g_{ii}^{(k)}−|g_{ij}^{(k)}| tan ψ^{(k)}} , \\ e^{ı ω_{k}} & = & \frac{z_{k}}{|z_{k}|}, i.e. ω_{k} = arg(z_{k}) , \\ tan φ^{(k)} & = & |z_{k}|, φ^{(k)} ∈ [0, \frac{π}{2} ]. \end{matrix} \}  (2.3)

Using convenient notation: c_{φ}^{(k)} = cos φ^{(k)}, s_{φ}^{(k)} = sin φ^{(k)}, c_{ψ}^{(k)} = cos ψ^{(k)},

s_{ψ}^{(k)} = sin ψ^{(k)}, the transformation formulas for the elements of G^{(k)} are given by the relation (2.4). The transformation G^{(k)} → G^{(k+1)}, defined by the relations (2.3) and (2.4) is referred to as Kogbetliantz transformation or Kogbetliantz step on the associated matrix G^{(k)}, with pivot pair (i, j). Shorter, we call it (associated) Kogbetliantz transformation or step. The algorithm (and the whole iterative process) on matrices G^{(k)} is called associated Kogbetliantz algorithm (process) or shorter Kogbetliantz algorithm (process) on G.

\begin{matrix}g_{li}^{(k+1)} = c_{φ}^{(k)}g_{li}^{(k)} − e^{ı ω_{k}}s_{φ}^{(k)}g_{lj}^{(k)} \\ g_{lj}^{(k+1)} = e^{−ı ω_{k}}s_{φ}^{(k)}g_{li}^{(k)} + c_{φ}^{(k)}g_{lj}^{(k)} \end{matrix} \} 1 ≤ l ≤ i − 1 ; \begin{matrix}g_{il}^{(k+1)} = c_{ψ}^{(k)}g_{il}^{(k)} − e^{−ı ϑ_{k}}s_{ψ}^{(k)}g_{lj}^{(k)} \\ g_{lj}^{(k+1)} = e^{ı ϑ_{k}}s_{ψ}^{(k)}g_{il}^{(k)} + c_{ψ}^{(k)}g_{lj}^{(k)} \end{matrix} \} i + 1 ≤ l ≤ j − 1 ; \begin{matrix}g_{il}^{(k+1)} = c_{φ}^{(k)}g_{il}^{(k)} − e^{ı ω_{k}}s_{φ}^{(k)}g_{jl}^{(k)} \\ g_{jl}^{(k+1)} = e^{−ı ω_{k}}s_{φ}^{(k)}g_{il}^{(k)} + c_{φ}^{(k)}g_{jl}^{(k)} \end{matrix} \} j + 1 ≤ l ≤ n ; \begin{matrix}g_{ii}^{(k+1)} = \frac{c_{φ}^{(k)}}{c_{ψ}^{(k)}} g_{ii}^{(k)} = \frac{|s_{ψ}^{(k)}|}{s_{φ}^{(k)}}g_{jj}^{(k)} = \frac{c_{ψ}^{(k)}}{c_{φ}^{(k)}} g_{ii}^{(k)} − \frac{s_{ψ}^{(k)}}{c_{φ}^{(k)}} |g_{ij}^{(k)}| \\ g_{jj}^{(k+1)} = \frac{c_{ψ}^{(k)}}{c_{φ}^{(k)}} g_{jj}^{(k)} = \frac{s_{φ}^{(k)}}{|s_{ψ}^{(k)}|}g_{ii}^{(k)} = \frac{c_{φ}^{(k)}}{c_{ψ}^{(k)}} g_{jj}^{(k)} + sgn (s_{ψ}^{(k)})\frac{s_{φ}^{(k)}}{c_{ψ}^{(k)}} |g_{ij}^{(k)}| \\ g_{ij}^{(k+1)} = 0 . \end{matrix}  (2.4)

Starting with upper-triangular B, the relation (2.4) holds for odd cycles only. For even cycles φ^{(k)}, ω_{k} and ψ^{(k)}, ϑ_{k} interchange their places. We can consider a method in which U^{(k)} and V^{(k)} interchange their places in subsequent cycles and thus ensure that (2.4) holds for every cycle.^{2} But we need not be concerned about that since we shall use estimates from [20] which are symmetric in the arguments φ^{(k)} and ψ^{(k)}. For this reason we do not consider the angle formulas when U^{(k)} is calculated first.

Let us now consider the scaled iterates of the matrices G^{(k)},

G_{S}^{(k)} = |diag(G^{(k)})|^{−1/2}G^{(k)}|diag(G^{(k)})|^{−1/2}, k ≥ 0.  (2.5)

The matrix G_{S}^{(k)} is associated with the content of the upper-triangular array when B_{S}^{(k)} is compactly stored in it. Next, we define

A^{(k)} = Ω(G_{S}^{(k)} + G_{S}^{(k)T}) = Ω(G_{S}^{(k)}) + Ω(G_{S}^{(k)})^{T} , k ≥ 0.  (2.6)

Since each B_{S}^{(k)} is essentially triangular, the relations (2.5) and (2.6) imply A^{(k)} = Ω(B_{S}^{(k)}) + Ω(B_{S}^{(k)})^{T} , k ≥ 0. Note that A^{(k)} = (a_{lm}^{(k)}) is a (complex) symmetric matrix with zero diagonal, whose off-diagonal elements are given by

a_{lm}^{(k)} = a_{ml}^{(k)} = \frac{g_{lm}^{(k)}}{\sqrt{g_{ll}^{(k)}g_{mm}^{(k)}}} = \frac{b_{lm}^{(k)} + b_{ml}^{(k)}}{\sqrt{b_{ll}^{(k)}b_{mm}^{(k)}}} , l ≠ m, k ≥ 0.  (2.7)

2The interchanges reflect only on singular vector updates.

The Frobenius norm of Ω(G_{S}^{(k)}) shall be denoted by α_{k},

α_{k} = ||Ω(G_{S}^{(k)})|| = ||Ω(B_{S}^{(k)})|| = \frac{\sqrt{2}}{2}||A^{(k)}||, k ≥ 0.  (2.8)

This paper estimates the quadratic reduction of α_{k} within one cycle. To this end, we can presume that k = 0.

3 Preliminary results

We recall several auxiliary results from [20] concerning nonorthogonal transformation G_{S}^{(k)} → G_{S}^{(k+1)} which do not depend on multiplicities of singular values. We assume that the initial matrix has positive diagonal elements.

Lemma 3.1 Suppose the sequence G^{(k)} = (g_{lm}^{(k)}), k ≥ 0 is obtained by the associated Kogbetliantz process, and suppose that G^{(0)} has positive diagonal elements. For a fixed k, let us denote G = (g_{lm}) = G^{(k)} and G̃ = (g̃_{lm}) = G^{(k+1)}. Thus, G̃ is obtained from G by a single Kogbetliantz transformation which annihilates the element g_{ij}. Let A = (a_{lm}) and Ã = (ã_{lm}) be defined by A = Ω(G_{S}) + Ω(G_{S})^{T} , Ã = Ω(G̃_{S}) + Ω(G̃_{S})^{T}, respectively.

(i) If |a_{ij}| < 1, then |ã_{li}|^{2} + |ã_{lj}|^{2} ≤ \frac{1 + |a_{ij}|}{1 − |a_{ij}|}(|a_{li}|^{2} + |a_{lj}|^{2}), l ≠ i, j .

(ii) If |a_{ij}| < 1 and ∥Ω(G̃_{S})∥ > ∥Ω(G_{S})∥, then

(a) ∥Ω(G̃_{S})∥^{2} − ∥Ω(G_{S})∥^{2} ≤ |a_{ij}| · \frac{2∥Ω(G_{S})∥^{2} − |a_{ij}|}{1 − |a_{ij}|} and

(b) |a_{ij}| ≤ 2∥Ω(G_{S})∥^{2} .

Proof. See the proof of [20, Lemma 1 and Lemma 2].

The next lemma bounds the growth of the off-norm of scaled matrices which can temporarily increase during the process.

Lemma 3.2 Let G^{(0)} = G, G^{(1)},…,G^{(N)} be obtained by applying N Kogbetliantz steps to G and let α_{k}, 0 ≤ k ≤ N be defined by the relation (2.8). If

α_{0} ≤ \frac{1}{10n}, n ≥ 3,

then

α_{k}^{2} ≤ c_{k}α_{0}^{2}, c_{k} = (1 + \frac{0.011}{n^{2}})^{k} < 1.006 0 ≤ k ≤ N .

Proof. See the proof of [20, Lemma 3].

The preceeding lemmas hold for both the column- and the row-cyclic strategies. The next lemma estimates the growth of the affected part of the rth column after s(s − 1)/2 Kogbetliantz steps under the column-cyclic strategy. Let

η_{sr}^{(k)} = [a_{1r}^{(k)}, …, a_{sr}^{(k)}]^{T} , 1 ≤ s < r ≤ n, 0 ≤ k ≤ N and
Q_{s} = 1 + 2 + … + (s − 1) .

Lemma 3.3 Let G^{(0)} = G, G^{(1)},…,G^{(N)} be obtained by applying N Kogbetliantz steps to G under the column-cyclic strategy.

If α_{0} ≤ \frac{1}{10n} then ∥η_{sr}^{(Q_{s})}∥^{2} ≤ K_{s}∥η_{sr}^{(0)}∥^{2}

where

K_{s} = [(1 + α)/(1 − α)]^{−(2s−3)} < 1.565 and α = \sqrt{1.006} α_{0} .  (3.1)

Proof. See the proof of [17, 18, 20, Lemma 4]. Let

a = max \{\frac{|g_{i(k)j(k)}^{(k)}|}{\sqrt{g_{i(k)i(k)}^{(k)}g_{j(k)j(k)}^{(k)}}} , 0 ≤ k < Q_{s} \} .

By Lemma 3.2, a ≤ \sqrt{1.006} α_{0} = α < 1. According to Lemma 3.1(i) instead of the term 1/(1 − a) in [18, Lemma 4], we now have (1 + a)/(1 − a). Following the lines in the proof of [18, Lemma 4], we obtain

||η_{sr}^{(Q_{s})}∥^{2} ≤ (\frac{1 + a}{1 − a})^{2s−3} ∥η_{sr}^{(0)}∥^{2} .  (3.2)

By the assumption on α_{0}, we have n · 2a/(1 − a) < 1. Therefore, we can use the inequality (1 + x)^{k} ≤ 1/(1 − kx). We obtain

(\frac{1 + a}{1 − a})^{2s−3} ≤ (1 + \frac{2a}{1 − a})^{2(n−1)} ≤ (1 − \frac{2(n − 1)a}{1 − a})^{−2}
= (1 + \frac{2(n − 1)a}{1 + a − 2na})^{2} ≤ (1 + \frac{2na}{1 − 2na})^{2} < 1.565 . (3.3)

4 Asymptotic assumptions

Here we define asymptotic assumptions for B = (b_{lm}) and prove some auxiliary results.

(B1) B ∈C^{n×n}, n ≥ 3 is upper triangular nonsingular matrix, satisfying

α_{0} ≤ \frac{1}{10} min\{\frac{1}{n}, θ\},

where α_{0} and θ are defined by the relations (2.8) and (1.6), respectively.

(B2) Diagonal elements of B are positive and satisfy

(∀l ∈ N_{q}) (∀m ∈ N_{r}) q < r ⇒ b_{ll} > b_{mm}

where the sets N_{r}, 1 ≤ r ≤ p are defined by the relation (1.2).

The singular values of B are given by the relation (1.1) and their relative gaps by the relations (1.5) and (1.6). Since B is nonsingular all relative gaps are smaller than one. It is clear that the assumptions (B1) and (B2) hold also for the matrix G (= B).

Lemma 4.1 Let G = G^{(0)} satisfy the asymptotic assumptions (B1) and (B2). Let the sequence G^{(k)} = (g_{lm}^{(k)}), k = 0, 1, … N be generated by the associated Kogbetliantz algorithm. Then for 0 ≤ k ≤ N holds

(i)  (1 − \frac{θ}{49})σ_{s_{r}} < g_{tt}^{(k)} < (1 + \frac{θ}{49})σ_{s_{r}}, t ∈ N_{r}, 1 ≤ r ≤ p ,
(ii)  rg (g_{ll}^{(k)}, g_{mm}^{(k)}) > \frac{24}{25}θ, l ∈ N_{q}, m ∈ N_{r}, q ≠ r ,
(iii)  | tan 2ψ^{(k)}| ≤ \frac{2|g_{ij}^{(k)}|}{g_{ii}^{(k)} − g_{jj}^{(k)}} , i(k) ∈ N_{q}, j(k) ∈ N_{r}, q < r,

(iv)  max\{| tan ψ^{(k)}|, | tan φ^{(k)}|\} ≤ \frac{|g_{ij}^{(k)}|}{g_{ii}^{(k)} − g_{jj}^{(k)}} ≤ \frac{25}{48} \frac{|a_{ij}^{(k)}|}{θ} ,

i(k) ∈ N_{q}, j(k) ∈ N_{r}, q < r.

where rg (·, ·) and a_{ij}^{(k)} are defined by the relations (1.4) and (2.7), respectively.

Proof. (i) We follow the proof of [18, Lemma 5(i)]. Note that (B1) and Lemma 3.2 ensure that condition (b) of Theorem 1.1 holds for all 0 ≤ k ≤ N. (B2) ensures that condition (a) of Theorem 1.1 holds for k = 0. We use the induction over k. By Theorem 1.1(i) we have for t ∈ N_{r}, 1 ≤ r ≤ p , and some k,

1 − \frac{2α_{k}^{2}}{θ_{r}} ≤ \frac{σ_{s_{r}}}{g_{tt}^{(k)}} ≤ 1 + \frac{2α_{k}^{2}}{θ_{r}} ,

or equivalently

σ_{s_{r}}(1 − \frac{\frac{2α_{k}^{2}}{θ_{r}}}{1 + \frac{2α_{k}^{2}}{θ_{r}}} ) ≤ g_{tt}^{(k)} ≤ σ_{s_{r}}(1 + \frac{\frac{2α_{k}^{2}}{θ_{r}}}{1 − \frac{2α_{k}^{2}}{θ_{r}}} ) .

Using Lemma 3.2 and the assumption (B1) we obtain

σ_{s_{r}}(1 − \frac{2 · 1.006}{100}θ) ≤ g_{tt}^{(k)} ≤ σ_{s_{r}}(1 + \frac{\frac{2·1.006}{100}}{1 − 2 · 1.006 · \frac{1}{10} · \frac{1}{30}} θ) ,

which yields g_{tt}^{(k)} ∈ D_{r}, t ∈ N_{r} where D_{τ} = {x : (1 − θ/49)σ_{s_{τ}} < x < (1 + θ/49)σ_{s_{τ}} }, 1 ≤ τ ≤ p. Using the definition (1.6) of θ we obtain for r ≠ q,

(1 − \frac{θ}{49}) max\{σ_{s_{q}}, σ_{s_{r}}\} − (1 + \frac{θ}{49}) min\{σ_{s_{q}}, σ_{s_{r}}\}

≥ |σ_{s_{q}} − σ_{s_{r}}| − \frac{θ}{49}(σ_{s_{q}} + σ_{s_{r}}) ≥ θ(σ_{s_{q}} + σ_{s_{r}}) − \frac{θ}{49}(σ_{s_{q}} + σ_{s_{r}})
= (σ_{s_{q}} + σ_{s_{r}}) · (1 − \frac{1}{49})θ = \frac{48}{49}θ(σ_{s_{q}} + σ_{s_{r}}) .

which yields

| x − y |≥ \frac{48}{49}θ(σ_{s_{q}} + σ_{s_{r}}) whenever x ∈ D_{q}, y ∈ D_{r}, q ≠ r .

Using the obtained relations, in the same way as in the proof of [18, Lemma 5(i)], we obtain for i = i(k) ∈ N_{q}, j = j(k) ∈ N_{r}, q < r,

|g_{ii}^{(k+1)} − g_{ii}^{(k)}| ≤ |g_{ij}^{(k)}| = |a_{ij}^{(k)}| · \sqrt{g_{ii}^{(k)}g_{jj}^{(k)}}
≤ \sqrt{1.006}α_{0} · (1 + \frac{θ}{49}) · \sqrt{σ_{s_{r}}σ_{s_{q}}}
≤ \sqrt{1.006} · \frac{θ}{10} · \frac{50}{49} · \frac{1}{2}(σ_{s_{r}} + σ_{s_{q}})
≤ \frac{3}{49} · (σ_{s_{r}} + σ_{s_{q}}) .

This yields implications g_{ii}^{(k)} ∈ D_{q} ⇒ g_{ii}^{(k+1)} ∈ D_{q} and g_{jj}^{(k)} ∈ D_{r} ⇒ g_{jj}^{(k+1)} ∈ D_{r}. The rest of the proof of (i) follows the proof of [18, Lemma 5(i)].

(ii) The proof follows the lines of the proof of [18, Lemma 5(ii)] and uses (i).

(iii) To simplify notation we omit the superscript k. By (ii), Lemma 3.2 and the assumption (B1), we have

\frac{\frac{|g_{ij}|^{2}}{g_{jj}}}{g_{ii} − g_{jj}} = \frac{\frac{|g_{ij}|^{2}}{g_{ii}g_{jj}}}{\frac{g_{ii}−g_{jj}}{g_{ii}+g_{jj}} · \frac{g_{ii}+g_{jj}}{g_{ii}}} = \frac{|a_{ij}|^{2}}{rg(g_{ii}, g_{jj}) · (1 + \frac{g_{jj}}{g_{ii}})} ≤ \frac{|a_{ij}|^{2}}{\frac{24}{25}θ · 1} < 1 ,

hence |g_{ij}|^{2}/(g_{ii} − g_{jj}) < g_{jj} and (see the relation (2.3))

tan 2ψ = \frac{2|g_{ij}|}{g_{jj} − g_{ii}} · \frac{g_{ii}}{\frac{g_{jj}^{2}−g_{ii}^{2}+|g_{ij}|^{2}}{g_{jj}−g_{ii}}} = \frac{2|g_{ij}|}{g_{jj} − g_{ii}} · \frac{g_{ii}}{g_{ii} + (g_{jj} − \frac{|g_{ij}|^{2}}{g_{ii}−g_{jj}})} ,

which yields the assertion. Note that the assumption (B2) yields g_{ii} > g_{jj}.

(iv) Using the assertion (iii), the first inequality for the angle ψ follows from | tan ψ| ≤ \frac{1}{2}| tan 2ψ|. For the angle φ we use (ii) and the assumption (B1),

\frac{|g_{ij} tan ψ|}{g_{ii} − g_{jj}} ≤ \frac{|g_{ij}|^{2}}{(g_{ii} − g_{jj})^{2}} = \frac{|a_{ij}|^{2}}{(\frac{g_{ii}−g_{jj}}{g_{ii}+g_{jj}})^{2} · \frac{(g_{ii}+g_{jj})^{2}}{g_{ii}g_{jj}}} ≤ \frac{|a_{ij}|^{2}}{(\frac{24}{25}θ)^{2} · 4} ≤ \frac{|a_{ij}|^{2}}{θ^{2}} < 1.

Hence

g_{jj} < g_{ii} − |g_{ij} tan ψ| ≤ g_{ii} − |g_{ij}| tan ψ

and therefore

| tan φ| ≤ |\frac{g_{jj}}{g_{ii} − |g_{ij}| tan ψ} | · | tan ψ| ≤ | tan ψ| ,

which is the first inequality for the angle φ. In the rest of the proof we use

|g_{ij}|/(g_{ii} − g_{jj}) ≤ |a_{ij}|/(2rg(g_{ii}, g_{jj})) and the assertion (ii).

5 The main Theorem

Here we state and prove the quadratic convergence result.

Theorem 5.1 Let B ∈C^{n×n} satisfy the asymptotic assumptions (B1) and (B2). Let the sequence B^{(0)} = B, B^{(1)}, …, B^{(N)} be generated by the serial Kogbetliantz method. Then

α_{N} ≤ \frac{13}{2} · \frac{α_{0}^{2}}{θ}

where α_{0}, α_{N} and θ are defined by the relations (2.8) and (1.6) respectively.

This result is quite analogous to the one for the Hermitian Jacobi method (see [18, 17, 16]). Note that the classical result from [11] (see also [2, 8, 21, 6])

∥Ω(B^{(N+1)})∥_{F} ≤ 4 \frac{Ω(B)∥_{F}^{2}}{δ} provided that Ω(B) < \frac{δ}{3},  (5.1)

where δ is minimum absolute gap in the set {−σ_{s_{1}}, …, −σ_{s_{p}}, σ_{s_{p}}, …, σ_{s_{1}}}, is either not valid or in praxis not applicable for tiny δ. In particular, if σ_{s_{p}} is tiny, then δ (≤ 2σ_{s_{p}}) is also tiny, while θ need not be tiny. Hence, the new result is especially applicable in cases when singular values cluster around zero.

Proof of Theorem 5.1

Because the row- and column-cyclic strategies are equivalent it suffices to prove the theorem for the column-cyclic strategy. The proof is based on the proof of [18, Theorem 6]. It uses the same technique: induction over the set {1, 2, …, p}.

Let A^{(k)} = (a_{ij}^{(k)}), k ≥ 0 be defined by the relation (2.6). For 1 ≤ t ≤ n we define the following submatrices of A^{(k)}, k ≥ 0 (we recall notation from [18]),

M_{t}^{(k)} = (a_{ij}^{(k)}), 1 ≤ i, j ≤ t, M_{t} = M_{t}^{(0)},
N_{t}^{(k)} = (a_{ij}^{(k)}), 1 ≤ i ≤ t, 1 ≤ j ≤ n, N_{t} = N_{t}^{(0)},
T_{t}^{(k)} = (τ_{ij}^{(k)}), 1 ≤ i, j ≤ t, T_{t} = T_{t}^{(0)}

where

τ_{ij}^{(k)} = \{ \begin{matrix}a_{ij}^{(k)}, & 1 ≤ i ≤ s_{r−1} < j ≤ min\{t, s_{r}\}, 2 ≤ r ≤ p \\ 0, & otherwise \end{matrix}

and certain blocks of A^{(k)}, k ≥ 0 (see the figure below):

F_{r}^{(k)} = (a_{ij}^{(k)}), 1 ≤ i ≤ s_{r−1} < j ≤ s_{r}, F_{r} = F_{r}^{(0)},
\overline{F} _{r}^{(k)} = (a_{ij}^{(k)}), s_{r−1} < j ≤ s_{r} < i ≤ n, \overline{F}_{r} = \overline{F} _{r}^{(0)},
G_{r}^{(k)} = (a_{ij}^{(k)}), 1 ≤ i ≤ s_{r}, max{i, s_{r−1}} < j ≤ s_{r}, G_{r} = G_{r}^{(0)} .

[figure: s_{r−1}; s_{r}; s_{r−1}; s_{r}; F_{r}; F_{r}^{T}; \overline{F} _{r}^{T}; \overline{F}_{r}; s_{r−1}; s_{r}; s_{r−1}; s_{r}; G_{r}; G_{r}^{T}]

Using induction with respect to r we shall prove the key relation

∥T_{s_{r}}^{(Q_{sr})}∥ ≤ C_{r} · \frac{∥N_{s_{r}}∥^{2}}{θ} , 1 ≤ r ≤ p,  (5.2)

where Q_{s_{r}} = 1 + 2 + … + (s_{r} − 1), and

C_{r} = 2.2 ∏_{i=1}^{r}(1 + 50\frac{∥G_{i}∥^{2}}{θ^{2}} )^{1/2}, 1 ≤ r ≤ p .  (5.3)

Since ∑_{i=1}^{r} ∥G_{i}∥^{2}/θ^{2} ≤ α_{0}^{2}/θ^{2}, the assumption (B1) implies

C_{r}^{2} ≤ \frac{2.2^{2}}{1 − 50 ∑_{i=1}^{r} \frac{∥G_{i}∥^{2}}{θ^{2}}} ≤ \frac{2.2^{2}}{1 − 50 · \frac{1}{100}} = 9.68, 1 ≤ r ≤ p .  (5.4)

The proof of the relation (5.2) is divided into three parts. We shall follow the proof of [18, Theorem 6] and briefly elaborate the proof of Theorem 5.1. \underline{The induction base (Part I)} We prove that the relation (5.2) holds for r = 1. It is obvious because T_{s_{1}}^{(Q_{s_{1}} )} is null-matrix. \underline{The induction hypothesis:} We assume that the relation (5.2) holds for some r − 1, 1 < r ≤ p, i.e.

∥T_{s_{r−1}}^{(I)}∥ ≤ C_{r−1} · \frac{∥N_{s_{r−1}}∥^{2}}{θ} , I := Q_{s_{r−1}} = 1 + 2 + … + (s_{r−1} − 1). (5.5)

The induction step (Parts II and III) Part II Annihilation of the elements in block F_{r} by columns. We shall prove

∥T_{s_{r−1}}^{(II)}∥ ≤ C_{r−1} · \frac{∥N_{s_{r}}∥^{2}}{θ} ,  (5.6)
∥F_{r}^{(II)}∥ ≤ 13.5\frac{∥N_{s_{r}}∥^{2}}{θ^{2}} ∥G_{r}∥ ,  (5.7)

where II := Q_{s_{r−1}} + s_{r−1}(s_{r} − s_{r−1}). Proof of the relation (5.6) We shall use the superscript w_{m} defined by

w_{t} = Q_{s_{r−1}} + (t − 1 − s_{r−1})s_{r−1}, s_{r−1} + 1 ≤ t ≤ s_{r} .  (5.8)

Lemma 5.2 Let B satisfy the assumptions of Theorem 5.1. Let w_{m} be defined by the relation (5.8). Then

(i) a_{km}^{(w_{m}+k)} = 0, 1 ≤ k ≤ s_{r−1},

(ii) |a_{lm}^{(w_{m}+k)}| ≤ \sqrt{1.006}(|a_{lm}^{(w_{m})}| + \frac{25}{24θ} ∑_{t=1}^{k} |a_{lt}^{(w_{m})}a_{tm}^{(w_{m}+t−1)}|),

1 ≤ k ≤ s_{r−1}, k < l ≤ m − 1,

(iii) |a_{lm}^{(w_{m}+k)}| ≤ \frac{1.045}{θ} ∑_{t=l+1}^{k} |a_{tm}^{(w_{m}+t−1)}|(|a_{lt}^{(w_{m})}|+ \frac{0.27}{θ}|a_{lm}^{(w_{m}+l−1)}a_{tm}^{(w_{m}+l−1)}|),

1 ≤ l < k ≤ s_{r−1},

(iv) |a_{lt}^{(w_{m}+s_{r−1})}| ≤ |a_{lt}^{(w_{m})}| + \frac{25}{48θ} |a_{lm}^{(w_{m}+l−1)}a_{tm}^{(w_{m}+l−1)}|, 1 ≤ l ≠ t ≤ s_{r−1},

(iv) |a_{lt}^{(w_{m}+s_{r−1})}| ≤ |a_{lt}^{(w_{m})}|+ \frac{25}{48θ} |a_{lm}^{(w_{m}+l−1)}a_{tm}^{(w_{m}+l−1)}|, 1 ≤ l ≤ s_{r−1} < t < m.

Proof. In the proof we shall omit the superscript w_{m}. The assertion (i) is obvious because a_{km}^{(k)} is annihilated pivot element. (ii) Using Lemma 4.1(ii), (iv) the proof becomes analogous to the proof of [18, Lemma 8(ii)]. The only difference arises in deriving the bound for g_{mm}^{(t−1)}/g_{mm}^{(k)}, 1 ≤ t ≤ k ≤ s_{r−1}. The bound will be obtained in the following way. Using Lemma 4.1(ii), we have

\frac{|g_{km}^{(k−1)}|^{2}}{(g_{kk}^{(k−1)})^{2} − (g_{mm}^{(k−1)})^{2}} = \frac{|a_{km}^{(k−1)}|^{2}}{\underbrace{\frac{(g_{kk}^{(k−1)})^{2} − (g_{mm}^{(k−1)})^{2}}{(g_{kk}^{(k−1)} + g_{mm}^{(k−1)})^{2}}}_{≥24θ/25>θ/2} · \underbrace{\frac{(g_{kk}^{(k−1)} + g_{mm}^{(k−1)})^{2}}{g_{kk}^{(k−1)}g_{mm}^{(k−1)}}}_{≥4}}

≤ \frac{1}{2} · \frac{|a_{km}^{(k−1)}|^{2}}{θ} < 1 .

From the above relation we have

(g_{mm}^{(k−1)})^{2} − (g_{kk}^{(k−1)})^{2} + |g_{km}^{(k−1)}|^{2} < 0 ,

which together with the relation (2.3) yields

| tan φ^{(k−1)}| ≤ \frac{g_{mm}^{(k−1)}}{g_{kk}^{(k−1)}} | tan ψ^{(k−1)}| < | tan ψ^{(k−1)}| ≤ 1  (5.9)

and

| tan ψ^{(k−1)}| ≤ \frac{g_{mm}^{(k−1)}}{g_{kk}^{(k−1)}} | tan φ^{(k−1)}| + \frac{|g_{km}^{(k−1)}|}{g_{kk}^{(k−1)}} ≤ | tan φ^{(k−1)}| + |a_{km}^{(k−1)}|. (5.10)

Using now the relations (2.4), (5.9), (5.10) and Lemma 4.1(iv), we have

\frac{g_{mm}^{(k−1)}}{g_{mm}^{(k)}} = \frac{cos φ^{(k−1)}}{cos ψ^{(k−1)}} = \sqrt{1 + \frac{tan^{2} ψ^{(k−1)} − tan^{2} φ^{(k−1)}}{1 + tan^{2} φ^{(k−1)}}}
< \sqrt{1 + tan^{2} ψ^{(k−1)} − tan^{2} φ^{(k−1)}}
≤ 1 + \frac{1}{2}| tan ψ^{(k−1)} + tan φ^{(k−1)}| · | tan ψ^{(k−1)} − tan φ^{(k−1)}|
≤ 1 + \frac{1}{2} · 2| tan ψ^{(k−1)}| · |a_{km}^{(k−1)}|
≤ 1 + \frac{25}{48} · \frac{|a_{km}^{(k−1)}|^{2}}{θ} ≤ 1 + \frac{25}{48} · \frac{α^{2}}{θ} , 1 ≤ k ≤ s_{r−1} .

Using the assumption (B1) and the relation (3.1), we obtain

\frac{g_{mm}^{(t−1)}}{g_{mm}^{(k)}} ≤ (1 + \frac{25α^{2}}{48θ})^{k−t+1} < (1 + \frac{25α^{2}}{48θ})^{n} < \frac{1}{1 − \frac{25}{48}nα · \frac{α}{θ}}

≤ \frac{1}{1 − \frac{25·1.006}{48·100}} < 1.006 , 1 ≤ t ≤ k ≤ s_{r−1}

which is the desired bound.

Using Lemma 4.1(ii) and Lemma 4.1(iv), and following the proof of [18, Lemma 8 (iii)], [18, Lemma 8 (iv)] and [18, Lemma 8 (v)], we obtain (iii), (iv) and (v), respectively.

Using Lemma 5.2, we can estimate the elements of the mth column prior and after annihilations. Let

η_{m}^{(k)} = [a_{1m}^{(k)}, …, a_{s_{r−1},m}^{(k)}], s_{r−1} + 1 ≤ m ≤ s_{r}, 0 ≤ k ≤ N , (5.11)
ρ_{m} = [a_{1m}, a_{2m}, …, a_{m−1,m}], s_{r−1} + 1 ≤ m ≤ s_{r} .  (5.12)

Lemma 5.3 Let B satisfy the assumptions of Theorem 5.1. Let w_{m}, η_{m}^{(w_{m})} and M_{s_{r−1}}^{(w_{m})} be defined by the relations (5.8), (5.11) and (5.12), respectively. Then

(i) ∑_{l=1}^{s_{r−1}} | a_{lm}^{(w_{m}+q_{l})} |^{2}≤ μ_{m}^{2}∥η_{m}^{(w_{m})}∥^{2}, 0 ≤ q_{l} < l,

(ii) ∑_{l=1}^{s_{r−1}} | a_{lm}^{(w_{m}+q_{l})} |^{2}≤ \frac{0.55μ_{m}^{2}}{θ^{2}} ∥η_{m}^{(w_{m})}∥^{2}, l ≤ q_{l} ≤ s_{r−1},

where

μ_{m} = \frac{\sqrt{1.006}}{1 − \frac{0.74}{θ}∥M_{s_{r−1}}^{(w_{m})}∥} .  (5.13)

Proof. Using the assertions (ii) and (iii) of Lemma 5.2, we can prove the assertions (i) and (ii) in the same way as it is done in [18, Lemma 9].

Lemma 5.4 Let B satisfy the assumptions of Theorem 5.1. Then for s_{r−1} + 1 ≤ m ≤ s_{r} hold

(i) ∥M_{s_{r−1}}^{(w_{m}+s_{r−1})}∥ ≤ ∥M_{s_{r−1}}^{(w_{m})}∥ + \frac{25\sqrt{2}μ_{m}^{2}∥η_{m}^{(w_{m})}∥^{2}}{48θ}

·[1 + \frac{\sqrt{0.55}}{θ}(∥M_{s_{r−1}}^{(w_{m})}∥ + \frac{0.39μ_{m}^{2}}{θ} ∥η_{m}^{(w_{m})}∥^{2})],

(ii) ∥T_{s_{r−1}}^{(w_{m}+s_{r−1})}∥ ≤ ∥T_{s_{r−1}}^{(w_{m})}∥ + \frac{25μ_{m}^{2}∥η_{m}^{(w_{m})}∥^{2}}{48θ}

·[1 + \frac{\sqrt{0.55}}{θ}(∥M_{s_{r−1}}^{(w_{m})}∥ + \frac{0.39μ_{m}^{2}}{θ} ∥η_{m}^{(w_{m})}∥^{2})],

where w_{m}, η_{m}^{(k)} and μ_{m} are defined by the relations (5.8), (5.11) and (5.13), respectively.

Proof. Using Lemma 5.2(iv) and Lemma 5.3, we just follow the lines in the proof of [18, Lemma 10].

Lemma 5.5 Let B satisfy the assumptions of Theorem 5.1. Let s_{r−1} + 1 ≤ q < m ≤ s_{r} and let w_{m}, η_{m}^{(k)} and μ_{m} be defined by the relations (5.8), (5.11) and (5.13). Then

(i) |a_{qm}^{(w_{m}+l)}| ≤ \sqrt{1.006}(|a_{qm}^{(w_{m})}| + \frac{25μ_{m}}{24θ} ∥η_{q}^{(w_{m})}∥ · ∥η_{m}^{(w_{m})}∥),

(ii) |a_{qm}^{(w_{m})}| = |a_{qm}^{(w_{q+1})}| ≤ \sqrt{1.006}(|a_{qm}^{(w_{q})}| + \frac{25μ_{q}}{24θ} ∥η_{q}^{(w_{q})}∥ · ∥η_{m}^{(w_{q})}∥).

Proof. Using Lemma 5.2(ii) and Lemma 5.3(i), the proof is analogous to the proof of [18, Lemma 11].

Lemma 5.6 Let B satisfy the assumptions of Theorem 5.1. If the hypothesis (5.5) holds, then

∥M_{s_{r−1}}^{(I)}∥ ≤ 9.57\frac{∥N_{s_{r−1}}∥^{2}}{θ} .

Proof. Following the proof of [18, Lemma 12], we apply Theorem 1.1(i) to bound ∑_{t=1}^{r−1} ∥A_{tt}^{(I)}∥^{2}. Since by Lemma 3.2, the matrices B^{(k)} and G^{(k)} are α-s.d.d. ( α = \sqrt{1.006} α_{0}) for all 0 ≤ k ≤ N, we have for 1 ≤ t ≤ p, 0 ≤ k ≤ N

\frac{2}{θ}∥A_{tt}^{(k)}∥ ≤ \frac{4}{θ^{2}} [∥A_{tt}^{(k)}∥^{2} + 2∥F_{t}^{(k)}∥^{2} + 2∥\overline{F} _{t}^{(k)}∥^{2}] ≤ \frac{8}{θ^{2}} α^{2} .  (5.14)

Using Lemma 3.2 and the assumption (B1), we obtain for 0 ≤ k ≤ N,

∥A_{tt}^{(k)}∥ ≤ \frac{1}{1 − \frac{2}{θ}∥A_{tt}^{(k)}∥} · \frac{4}{θ}[∥F_{t}^{(k)}∥^{2} + ∥\overline{F} _{t}^{(k)}∥^{2}] ≤ 1.088 · \frac{4}{θ}[∥F_{t}^{(k)}∥^{2} + ∥\overline{F} _{t}^{(k)}∥^{2}].

(5.15)

Using the induction hypothesis (the relation (5.5)), Lemma 3.3 and the relation

\frac{∥N_{s_{r−1}}∥^{2}}{θ^{2}} ≤ \frac{2α_{0}^{2}}{θ^{2}} ,  (5.16)

we obtain (cf. the proof of [18, Lemma 12])

∥M_{s_{r−1}}∥^{2} ≤ [2C_{r−1}^{2} + 1.088^{2} · 16(4\frac{α_{0}^{2}}{θ^{2}} C_{r−1}^{2} + 1.565)^{2}] · \frac{∥N_{s_{r−1}}∥^{4}}{θ^{2}}

which, together with the relation (5.4), implies the assertion.

Lemma 5.7 Let B satisfy the assumptions of Theorem 5.1. Let w_{m} and μ_{m} be defined by the relations (5.8) and (5.13), respectively. If the hypothesis (5.5) holds, then for s_{r−1} + 1 ≤ m ≤ s_{r}, we have (i) μ_{m} ≤ μ = 1.169,

(ii) ∥M_{s_{r−1}}^{(w_{m+1})}∥ ≤ 9.57\frac{∥N_{m}∥^{2}}{θ},

(iii) ∥T_{s_{r−1}}^{(w_{m+1})}∥ ≤ C_{r−1} \frac{∥N_{m}∥^{2}}{θ}.

Proof. See the proof of [18, Lemma 13]. Induction with respect to m, s_{r−1}+1 ≤ m ≤ s_{r} is used.

(i) The relations (5.13), (5.16), Lemma 5.6 and the assumption (B1) are used. (ii) The assertion (i), the lemmas 5.4(i), 3.3, 5.6, the assumption (B1) and the relation (5.16) are used. (iii) Lemma 5.4(ii) and the relations (5.5) and (5.4) are used.

The assertion (iii) of Lemma 5.7 implies (by taking m = s_{r}) the relation (5.6).

Proof of the relation (5.7)

Lemma 5.8 Let B satisfy the assumptions of Theorem 5.1. Let s_{r−1} + 1 ≤ q ≤ s_{r} and let w_{m}, η_{m}^{(k)} and ρ_{m} be defined by the relations (5.8), (5.11) and (5.12), respectively. If the hypothesis (5.5) holds, then

∥η_{q}^{(w_{sr +1})}∥ ≤ 1.012∥η_{q}^{(w_{q+1})}∥ + \frac{0.774}{θ} ∑_{m=q+1}^{s_{r}} |a_{qm}^{(w_{q+1})}| · ∥ρ_{m}∥ .

Proof. Following the lines in the proof of [18, Lemma 14] and using Lemma 5.2 (v), Lemma 5.3(i), Lemma 5.5(i), Lemma 5.7(i) and Lemma 3.3, we obtain

∥η_{q}^{(w_{m+1})}∥ ≤ z_{m}∥η_{q}^{(w_{m})}∥ + z|a_{qm}^{(w_{q+1})}| · ∥ρ_{m}∥ , q + 1 ≤ m ≤ s_{r},

where

z_{m} = 1 + \frac{25 · 25 · \sqrt{1.006} · μ^{2} · 1.565}{48 · 24} · \frac{∥ρ_{m}∥^{2}}{θ^{2}} ≤ 1 + 1.164\frac{∥ρ_{m}∥^{2}}{θ^{2}} ,

z = \frac{25μ\sqrt{1.006 · 1.565}}{48θ} ≤ \frac{0.764}{θ} .

The rest of the proof is analogous to the proof of [18, Lemma 14] with the only exception that ∑_{m=1}^{n} ∥ρ_{m}∥^{2} is now bounded by α_{0}^{2}/θ^{2} (instead of α_{0}^{2}/(2γ^{2})).

Lemma 5.9 Let B satisfy the assumptions of Theorem 5.1. If the hypothesis (5.5) holds, then the relation (5.7) holds too.

Proof. See the proof of [18, Lemma 15]. Using the relation (5.14), we have

∥A_{rr}∥ ≤ \frac{4}{θ}∥N_{s_{r}}∥^{2} .  (5.17)

Now, using Lemma 5.3(ii), Lemma 5.7(i), (ii) and Lemma 3.3, we obtain

∥η_{q}^{(w_{q+1})}∥^{2} ≤ \frac{0.55 · 1.169^{2}}{θ^{2}} · 1.565∥ρ_{q}∥^{2}(9.57\frac{∥N_{q}∥^{2}}{θ})^{2}
≤ 107.73 · \frac{∥ρ_{q}∥^{2}}{θ^{2}} · \frac{∥N_{s_{r}}∥^{4}}{θ^{2}} , s_{r−1} + 1 ≤ q ≤ s_{r} .

Using Lemma 5.5(ii), Lemma 5.7(i) and Lemma 3.3, we have

∑_{q=s_{r−1}+1}^{s_{r}} ∑_{m=q+1}^{s_{r}} | a_{qm}^{(w_{q+1})} |^{2}≤ \sqrt{1.006}(\frac{∥A_{rr}∥}{\sqrt{2}} + \frac{25μ · 1.565}{48θ}∥G_{r}∥^{2})^{2} .

Combining the obtained expressions with Lemma 5.8 and the relation (5.17) as in the proof of [18, Lemma 15], we obtain the desired assertion.

Part III Annihilation of the elements in block A_{rr} by columns.

We shall prove

T_{s_{r−1}}^{(III)} = T_{s_{r−1}}^{(II)} ,  (5.18)
∥F_{r}^{(III)}∥ ≤ 1.07∥F_{r}^{(II)}∥ ,  (5.19)

where III := Q_{s_{r}} = 1 + 2 + … + (s_{r} − 1).

The relation (5.18) is obvios because the matrix T_{s_{r−1}}^{(II)} does not change at this stage. For s_{r−1} + 2 ≤ m ≤ s_{r} we define index

v_{m} = Q_{s_{r−1}} + (s_{r} − s_{r−1})s_{r−1} + (m − s_{r−1} − 1)(m − s_{r−1} − 2)/2 . (5.20)

The index v_{m} corresponds to the stage just before the annihilations in the mth column within A_{rr} take place.

Lemma 5.10 Let B satisfy the assumptions of Theorem 5.1. Then

∥F_{r}^{(v_{sr +1})}∥ ≤ 1.07 · ∥F_{r}^{(v_{s_{r−1}+2})}∥,

where v_{m} is defined by the relation (5.20).

Proof. See the proof of [18, Lemma 16]. Using the relation (5.14), we have

a ≤ max_{k,m} \frac{∥A_{rr}^{(v_{m}+k)}∥}{\sqrt{2}} ≤ \frac{1}{\sqrt{2}} · \frac{4}{θ} · α^{2} .

Using the relations (3.2) and (3.3), we obtain

(\frac{1 + a}{1 − a})^{2s−3} ≤ ( 1 + \frac{\sqrt{2} · 4 · nα · α/θ}{1 − \sqrt{2} · 4 · nα · α/θ} )^{2}

< (1 + \frac{\sqrt{2} · 4 · 1.006 · 0.01}{1 − \sqrt{2} · 4 · 1.006 · 0.01})^{2} < 1.07^{2} ,

which together with (3.1) and the assumption (B1) proves the lemma. Note that Lemma 5.10 proves the relation (5.19). This completes Part III. Using now the relations (5.18), (5.19), (5.6), (5.7) and (5.3), the induction step is quite analogous to that in [18]. Completing the proof of Theorem 5.1 Using the relation (5.15) we have

∑_{t=1}^{p} ∥A_{tt}^{(N)}∥^{2} ≤ 1.088^{2} · \frac{16}{θ^{2}} ∑_{t=1}^{p}[∥F_{t}^{(N)}∥^{2} + ∥\overline{F} _{t}^{(N)}∥^{2}]^{2}
≤ 1.088^{2} · \frac{16}{θ^{2}} · ∥T_{s_{p}}^{(N)}∥^{2} ∑_{t=1}^{p}[∥F_{t}^{(N)}∥^{2} + ∥\overline{F} _{t}^{(N)}∥^{2}]^{2}
≤ 1.088^{2} · \frac{16}{θ^{2}} · 2∥T_{s_{p}}^{(N)}∥^{4} < 37.88 · \frac{∥T_{s_{p}}^{(N)}∥^{4}}{θ^{2}} .

Finally, using the relation (5.4), one obtains

α_{N}^{2} = ∥T_{s_{p}}^{(N)}∥^{2} + \frac{1}{2} ∑_{t=1}^{p} ∥A_{tt}^{(N)}∥^{2} ≤ ∥T_{s_{p}}^{(N)}∥^{2} + \frac{1}{2} · 37.88 · \frac{∥T_{s_{p}}^{(N)}∥^{4}}{θ^{2}}

≤ ∥T_{s_{p}}^{(N)}∥^{2}(1 + 18.94 · \frac{∥T_{s_{p}}^{(N)}∥^{2}}{θ^{2}} ) ≤ C_{p}^{2} \frac{∥N_{s_{p}}∥^{4}}{θ^{2}} (1 + 18.94 · C_{p}^{2} · \frac{∥N_{s_{p}}∥^{4}}{θ^{4}} )
≤ 9.68 · \frac{(2α_{0}^{2})^{2}}{θ^{2}} (1 + 18.94 · 9.68 · \frac{(2α_{0}^{2})^{2}}{θ^{4}} ) ≤ 38.72(1 + 183.34 · 4 · 0.1^{4}) · \frac{α_{0}^{4}}{θ^{2}}
< (\frac{13}{2} · \frac{α_{0}^{2}}{θ} )^{2} ,

which completes the proof of Theorem 5.1.

6 Some Numerical Examples

We have performed some numerical tests. Test matrices have been generated in several ways and here we describe the simplest one. For given n, the uppertriangular matrix R(n) = (r_{ij}^{(n)}) is generated by “wrapping” the sequence of natural numbers: 1, 2, 3,… onto the superdiagonals of R(n), starting with the last one. Thus, r_{1n}^{(n)} = 1, r_{1,n−1}^{(n)} = 2, r_{2n}^{(n)} = 3,…,r_{11}^{(n)} = M_{n−1} + 1,…, r_{nn}^{(n)} = M_{n}, where M_{n} = n(n + 1)/2. Then R(n) is scaled as T (n) = D(n)R(n), T (n) = R(n)D(n) or T (n) = [D(n)]^{\frac{1}{2}} R(n)[D(n)]^{\frac{1}{2}} where diagonal matrix D(n) is defined by the vector d(n). Then B = T (n) is our initial matrix. The role of the transformation R(n) ↦ T (n) is to make T (n) scaled almost diagonal which means that the minimum absolute gap δ will be in most cases much smaller than the minimum relative gap θ. The vector d(n) is generated using four integer parameters: M, I1, I2, I3. The following FORTRAN 77 code has been used in obtaining d(n)

  DO I=1,M
   D(I)=1D1**(DBLE(I1)+DBLE((I2-I1)*(I-1))/DBLE(M-1))
  ENDDO
  DO I=M+1,N
   D(I)=1D1**(DBLE(I2)+DBLE((I3-I2)*(I-M))/DBLE(N-M))
  ENDDO

The parameters I1, I2 and I3 are integers between -20 and 8 and and M is between 1 and N, N is the matrix order. The parameters I1, I2, I3 have been varying by three nested loops, M has been set to N/2 and N has been varying by several values between 30 and 300.

The computation was done using double precision arithmetic. Visual Fortran Standard Edition 5.0D compiler of Digital corp. has been used.

We have considered the behavior of ∥Ω(B^{(rN)})∥_{F} and ∥Ω(B_{S}^{(rN)})∥_{F}, r = 1, 2, … and how well these measures qualify as predictors on the number of cycles needed to reach the stopping criterion. As stopping criterion, we have chosen the one of Rutishauser [23]. However, instead of using sophisticated Rutishauser’s transformation formulas, we have used the LAPACK routine DLASV2 for computing the transformation parameters and the BLAS 1 routine DROT for performing transformation.

In all performed diagonalizations, the quantity ∥Ω(B_{S}^{(rN)})∥_{F} has decreased as expected, smoothly quadratically, while ∥Ω(B^{(rN)})∥_{F} has often behaved unexpectedly.

Say, for N= 30, I1= 4, I2=-10, I3=-10, we have obtained d(30) = (10^{4}, 10^{3}, 10^{2}, 10, 1, 10^{−1}, 10^{−2}, 10^{−3}, 10^{−4}, 10^{−5}, 10^{−6}, 10^{−7}, 10^{−8}, …, 10^{−8}) and the starting upper-triangular matrix B = R(30)D(30) whose three largest and smallest diagonals and singular values are displayed in Table 1.

The absolute and relative gaps have been computed a posterior: δ ≈ 4.041 · 10^{−10}, θ ≈ 8.464 · 10^{−3}. Since θ >> δ, one can expect that the scaled iterates behave more regularly than the ordinary iterates. Table 2 below shows that straightforwardly.

We see that due to tiny δ, the quantity ∥Ω(B^{(rN)})∥_{F} behaves strange: strong reduction in the first cycle and no further quadratic reduction in the subsequent cycles.

\begin{aligned}i & b_{ii} & σ_{s_{i}} \\ 1 & 0.436E + 07 & ≈ 0.4379E+07 \\ 2 & 0.437E + 06 & ≈ 0.4370E+06 \\ 3 & 0.438E + 05 & ≈ 0.4380E+05 \\ 28 & 0.463E − 07 & ≈ 0.2463E-07 \\ 29 & 0.464E − 07 & ≈ 0.2407E-07 \\ 30 & 0.465E − 07 & ≈ 0.2367E-07 \end{aligned}

Table 1

\begin{matrix}r & ∥Ω(B^{(rN)})∥_{F} & ∥Ω(B_{S}^{(rN)})∥_{F} \\ 1 & 0.358E + 03 & 0.224E + 01 \\ 2 & 0.781E − 07 & 0.499E + 00 \\ 3 & 0.313E − 08 & 0.533E − 01 \\ 4 & 0.364E − 10 & 0.772E − 03 \\ 5 & 0.682E − 14 & 0.151E − 06 \\ 6 & 0.234E − 21 & 0.488E − 14 \\ 7 & 0.173E − 36 & 0.322E − 29 \end{matrix} Table 2

On the contrary, due to large θ, ∥Ω(B_{S}^{(rN)})∥_{F} exhibits sound quadratic reduction after second cycle. Not only that ∥Ω(B_{S}^{(rN)})∥_{F} reduces quadratically, but it gives a lower bound on the number of sure digits in diagonal elements as approximations of singular values. Even better bound for that is the quantity 2∥Ω(B_{S}^{(rN)})∥_{F}^{2} /θ (see Theorem 1.1(i)). Note that in cases like here, when ∥Ω(B_{S}^{(rN)})∥_{F} < θ/3, there are simple and cheap ways to compute the needed lower bound for all θ_{i}, 1 ≤ i ≤ p and hence for θ (see [19]).

The clue for strange behavior of ∥Ω(B^{(rN)})∥_{F} lies in the fact that the classical result (5.1) ensures the quadratic convergence of ∥Ω(B^{(rN)})∥_{F} /δ and not of ∥Ω(B^{(rN)})∥_{F} alone, provided that it is smaller than one.

However, if δ = 2σ_{s_{p}} and σ_{s_{p}} is tiny, we just have to blame the result (5.1) or the definition of δ which depends on the minimum singular value. In such a case, the classical result is misleading and useless. On the contrary, the new result solves the problem. Since θ does not depend on σ_{s_{p}} and since θ can be large even if singular values cluster around zero, the quadratic convergence ∥Ω(B_{S}^{(rN)})∥_{F} can be noticed and used even in early stages of the process.

We end the paper by another four examples of matrices. For each diagonalization process we display a table with relevant data.

    N= 30, I1=8, I2=-10 I3=-10

\begin{matrix}r & ∥Ω(B^{(rN)})∥_{F} & ∥Ω(B_{S}^{(rN)})∥_{F} \\ 1 & 0.108E + 09 & 0.224E + 01 \\ 2 & 0.476E + 02 & 0.498E + 00 \\ 3 & 0.313E − 08 & 0.533E − 01 \\ 4 & 0.364E − 10 & 0.772E − 03 \\ 5 & 0.682E − 14 & 0.151E − 06 \\ 6 & 0.234E − 21 & 0.488E − 14 \\ 7 & 0.174E − 36 & 0.322E − 29 \end{matrix}

                                           N= 150, I1=8, I2=-20 I3=-20

\begin{matrix}r & ∥Ω(B^{(rN)})∥_{F} & ∥Ω(B_{S}^{(rN)})∥_{F} \\ 1 & 0.364E + 13 & 0.910E + 01 \\ 2 & 0.877E + 10 & 0.236E + 01 \\ 3 & 0.149E + 07 & 0.530E + 00 \\ 4 & 0.123E − 01 & 0.672E − 01 \\ 5 & 0.347E − 18 & 0.161E − 02 \\ 6 & 0.221E − 21 & 0.100E − 05 \\ 7 & 0.672E − 28 & 0.342E − 12 \\ 8 & 0.673E − 41 & 0.475E − 25 \end{matrix}

δ = 4.04 · 10^{−10}, θ = 8.46 · 10^{−3}

Table 3: B= D^{\frac{1}{2}} R(N)D^{\frac{1}{2}}

δ = 8.21 · 10^{−20}, θ = 7.21 · 10^{−4}

Table 4: B = DR(N)

N= 150, I1=2, I2=-20 I3=-15

\begin{matrix}r & ∥Ω(B^{(rN)})∥_{F} & ∥Ω(B_{S}^{(rN)})∥_{F} \\ 1 & 0.580E + 06 & 0.288E + 03 \\ 2 & 0.841E + 04 & 0.189E + 01 \\ 3 & 0.360E + 00 & 0.556E + 00 \\ 4 & 0.150E − 09 & 0.153E + 00 \\ 5 & 0.179E − 13 & 0.187E − 01 \\ 6 & 0.397E − 15 & 0.412E − 03 \\ 7 & 0.439E − 19 & 0.420E − 07 \\ 8 & 0.168E − 27 & 0.140E − 15 \end{matrix}

δ = 8.35683 · 10^{−19}, θ = 5.27 · 10^{−3} Table 5: B = R(N)D

                                       N= 300, I1=6, I2=-20 I3=-15

\begin{matrix}r & ∥Ω(B^{(rN)})∥_{F} & ∥Ω(B_{S}^{(rN)})∥_{F} \\ 1 & 0.456E + 12 & 0.110E + 02 \\ 2 & 0.656E + 10 & 0.188E + 01 \\ 3 & 0.151E + 08 & 0.589E + 00 \\ 4 & 0.302E + 04 & 0.139E + 00 \\ 5 & 0.464E − 04 & 0.196E − 01 \\ 6 & 0.267E − 15 & 0.862E − 03 \\ 7 & 0.253E − 18 & 0.853E − 06 \\ 8 & 0.123E − 23 & 0.406E − 11 \\ 9 & 0.365E − 35 & 0.118E − 22 \end{matrix}

δ = 8.75 · 10^{−19}, θ = 1.63 · 10^{−3}

Table 6: B = DR(N) .

References

[1] Barlow J., Demmel J. W. (1990) Computing Accurate Eigensystems of Scaled Diagonally Dominant Matrices , SIAM J. Num. Anal. 27, 762-791.

[2] Charlier J. P., Van Dooren P. (1987) On Kogbetliantz’s SVD Algorithm in the Presence of Clusters, Linear Algebra Appl. 95, 135-160.

[3] Charlier J. P., Vanbegin M., Van Dooren P. (1988) On Efficient Implementation of Kogbetliantz’s Algorithm for Computing the Singular Value Decomposition, Numer Math. 52, 279-300.

[4] Demmel J. W., Veselić K. (1992) Jacobi’s Method is more Accurate then QR, SIAM J. Matrix Anal. Appl. 13, 1204-1245.

[5] Drmač Z. (1999): A posteriori computation of the singular vectors in a preconditioned Jacobi SVD algorithm, IMA J. Numer. Anal. 19, 191-213.

[6] Fernando K. V. (1989) Linear Convergence of the Row Cyclic Jacobi and Kogbetliantz Methods. Numer. Math. 56, 73–94.

[7] Forsythe G. E., Henrici P. (1960) The cyclic Jacobi method for computing the principal values of a complex matrix, Trans. Amer. Math. Soc. 94, 1–23.

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

[9] Hari V. and Veselić K. (1987) On Jacobi Methods for Singular Value Decompositions, SIAM J. Sci. Stat. Comput. Vol. 8, No. 5, 741–754.

[10] Hari V. (1989) On the Quadratic Convergence Bounds of the Serial Singular Value Decomposition Jacobi Methods for Triangular Matrices, SIAM J. Sci. Stat. Comput. 10, 1076-1096.

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

[12] Hari V., Drmač Z. (1995) On Scaled Almost Diagonal Hermitian Matrix Pairs, SIAM J. Matrix Anal. Appl. Vol. 18. No. 4 (1997) 1000-1012.

[13] Hari V. (1999) Parallel SVD methods for matrices in butterfly form, Technical report, University of Zagreb, (Croatian language).

[14] Kogbetliantz E. (1954) Diagonalization of General Complex Matrices as a New Method for Solution of Linear Equations, Proc. Imtern. Congr. Math. Amsterdam 2, 356–357.

[15] Kogbetliantz E. (1955) Solutions of Linear Equations by Diagonalization of Coeficient Matrices, Quart. Appl. Math. 13, 123–132.

[16] Matejaš J. (1996) Quadratic Convergence of Scaled Iterates by Diagonalization Methods, University of Zagreb, Ph. D. Thesis (Croatian language).

[17] Matejaš J. (1999) Quadratic Convergence of Scaled Matrices in Jacobi Method, Numer. Math. 87, 171-199.

[18] Matejaš J. (2000) Convergence of Scaled Iterates in Jacobi Method, (to appear in Linear Algebra and Its Applications).

[19] Matejaš J., Hari V. (1998) Scaled Almost Diagonal Matrices with Multiple Singular Values, ZAMM 78, 121-131.

[20] Matejaš J., Hari V. (1999) Scaled Iterates by Kogbetliantz Method, Proceedings of the Conference on Applied Mathematics and Computation, Dubrovnik 13-18. Sept. 1999, pp. 1-20.

[21] Paige C., Van Dooren, P. (1986) On the Quadratic Convergence of Kogbetliantz’s Algorithm for Computing the Singular Value Decomposition, Linear Algebra and Its Applications 77, 301–313.

[22] Voevodin V.V. (1966) Cislennye metody linejnoj algebry, Nauka, Moscow (in Russian).

[23] Wilkinson J.H. and Reinsch C. (1971) Handbook of Automatic Computation II, Linear Algebra, Springer, Berlin Heidelberg New York.