Contemporary Mathematics

Approximate displacement rank and applications

Dario Andrea Bini and Beatrice Meini

Abstract. We introduce the concepts of orthogonal displacement representation and of ϵ-displacement rank of matrix, prove some properties and show some applications to Toeplitz-like matrix inversion by means of Newton’s iteration, and to cyclic reduction for the solution of block banded block Toeplitz systems.

1. Introduction and definitions

The concepts of displacement rank and displacement operator, introduced by Kailath et al. [18, 17], and elaborated in many other papers (see, for instance, [15, 8, 22] and the references cited therein), provide powerful tools for dealing with Toeplitz matrices, allow the generalization of the latter structure, and are the basis for the design of inversion formulae and of fast inversion algorithm for structured matrices. We recall that a Toeplitz matrix has entries that are constant along the diagonals.

Different displacement operators have been introduced and analyzed in the literature [15, 8, 22, 19, 20, 21, 7], most of them are strongly related to representation formulae for Toeplitz-like matrices and for their inverses. Operators of the kind Δ_{S}(X) = XS − SX defined for X ∈ ℂ^{n×n}, where S ∈ ℂ^{n×n} is a given matrix, are clearly singular and their kernel is the centralizer of S. Any matrix X can be decomposed as sum of an element of the kernel and one of the image of Δ_{S}. Suitable choices of S provide different operators well suited for dealing with matrices of the Toeplitz type.

For instance, if S = Z, where Z = (z_{i,j}) is the n × n lower triangular Toeplitz matrix such that z_{i+1,i} = 1, z_{i,j} = 0 for i ≠ j + 1, then Δ_{Z}(X) is zero for any lower triangular Toeplitz matrix X, moreover, rankΔ_{Z}(A) ≤ 2 for any Toeplitz matrix

\begin{matrix}c \\ ◯ \end{matrix}0000 (copyright holder)

1991 Mathematics Subject Classification. Primary ; Secondary .

Supported in part by NSF Grant #000000.

A. The following well known results hold (see for instance [8])

Δ_{Z}(AB) = Δ_{Z}(A)B + AΔ_{Z}(B),
Δ_{Z}(A^{−1}) = −A^{−1}Δ_{Z}(A)A^{−1},
A = L(Ae_{1}) + ∑_{i=1}^{k} L(u_{i})L(Zv_{i})^{T} iff AZ − ZA = ∑_{i=1}^{k} u_{i}v_{i}^{T} ,
A^{−1} = L(A^{−1}e_{1}) − ∑_{i=1}^{k} L(A^{−1}u_{i})L(ZA^{−T} v_{i})^{T} iff
AZ − ZA = ∑_{i=1}^{k} u_{i}v_{i}^{T} , det A ≠ 0,

(1.1)

where L(x) = ∑_{i=1}^{n} x_{i}Z^{i−1} is the lower triangular Toeplitz matrix having x as first column, and e_{1} is the first column of the identity matrix. A matrix that can be written as sum of k products of lower and upper triangular Toeplitz matrices is called Toeplitz-like if k is independent of n. The representation of A given in terms of triangular Toeplitz matrices is called displacement representation.

We call displacement rank of A and denote it as drk(A), the rank of Δ_{Z}(A). In this way, a matrix with displacement rank k as well as its inverse (if nonsingular) can be written as sum of at most k + 1 products of lower and upper triangular Toeplitz matrices.

In the original formulation of [17, 18] the nonsingular displacement operator X − ZXZ^{T} has been used. This operator can be easily related to Δ_{Z}(A).

Similar operators and similar relations can be constructed by replacing Z with suitable matrices. It is worth citing the case S = H = Z + Z^{T}, that provides an operator suitable for Toeplitz plus Hankel matrices [7], or the choice S = C = Z +e_{1}e_{n}^{T} (unit circulant) that is still convenient for manipulating Toeplitz matrices. Here e_{n} denotes the last column of the n × n identity matrix.

Although the displacement rank is kept under inversion, it is not maintained under addition and multiplication. This is an unpleasant but unavoidable situation. Consider for instance the celebrated Newton iteration

\begin{matrix}X_{i+1} = 2X_{i} − X_{i}AX_{i}, i = 0, 1, 2, … \\ X_{0} = A^{T} θ, θ = 1/||AA^{T} || \end{matrix} (1.2)

such that lim_{i} X_{i} = A^{−1}, for a nonsingular matrix A. Here and hereafter || · || denotes any operator norm, and, if not specified, || · || is the Euclidean norm. If A is Toeplitz, then the displacement rank of X_{0} is at most 2, but the displacement rank of X_{i} grows exponentially with i until the value of n is reached. This fact nullifies the computational advantages of the displacement representation so that Newton’s iteration is not used in practice for Toeplitz matrix inversion for the large computational cost. In this case the displacement rank is not a useful tool for the computational treatment of the problem. An attempt to use Newton’s iteration for Toeplitz inversion is done in [23].

Another drawback of displacement representation (and inversion) formulae is related to their numerical stability. Observe that, since the decomposition of a matrix as sum of k matrices of rank 1 is not unique, then also the displacement representation formula of (1.1) is not unique. Moreover, if some vectors u_{i} and v_{i} are very large in norm with respect to the size of the entries of A, using the displacement representation formula of (1.1) for computing a matrix-vector product surely leads to numerical cancellation problems with loss of accuracy in the computed result. Issues related to the stability of such formulae are investigated in [14]. Indeed, a formula constructed with vectors u_{i}, v_{i} that have norms comparable with the norm of A is more desirable.

In this paper we introduce the concepts of orthogonal displacement representation and of approximate (ϵ-)displacement rank of a matrix that somehow weakens the concept of displacement rank and provides a more suitable tool for the numerical treatment of Toeplitz matrices. We show that certain numerical computations, like Newton’s iteration applied to Toeplitz-like matrices and cyclic reduction for block banded block Toeplitz matrices with Toeplitz blocks [4], can be performed at a very low computational cost, due to the properties of the ϵ-displacement rank. Moreover, by using the orthogonal displacement representation we obtain formulae that avoid the drawback of numerical cancellation.

For notational simplicity we write Δ(A) in place of Δ_{Z}(A)

Definition 1.1. For a given ϵ > 0 define the absolute ϵ-displacement rank of a matrix A as

drk_{ϵ}(A) = min_{||E||≤ϵ} rank(Δ(A) + E),

and the relative ϵ-displacement rank as

rdrk_{ϵ}(A) = min_{||E||≤ϵ||A||} rank(Δ(A) + E),

where || · || denotes the Euclidean norm.

If r = drk_{ϵ}(A) (r = rdrk_{ϵ}(A), respectively), and if Δ(A) = ∑_{i=1}^{r} u_{i}v_{i}^{T} + E, with ||E|| ≤ ϵ (||E|| ≤ ϵ||A||, respectively), from the above definition and from (1.1), we find that the matrix

A_{ϵ} = L(Ae_{1}) + ∑_{i=1}^{r} L(u_{i})L(Zv_{i})^{T}

provides an approximation to A given in terms of the approximate displacement representation. We refer to A_{ϵ} as to an ϵ-displacement (relative ϵ-displacement) representation of A, and we call the vectors u_{i}, v_{i}, i = 1, …, r, ϵ-displacement (relative ϵ-displacement) generators for A.

2. Main properties of ϵ-displacement rank

Let

Δ(A) = UΣV^{T} = ∑_{i=1}^{k} σ_{i}u_{i}v_{i}^{T} , σ_{1} ≥ σ_{2} ≥ … ≥ σ_{k} > 0 (2.1)

be the singular values decomposition (SVD) of Δ(A). We define orthogonal displacement representation (odr) of A (with respect to Δ) the expression

A = L(Ae_{1}) + ∑_{i=1}^{k} σ_{i}L(u_{i})L(Zv_{i})^{T} , (2.2)

and orthogonal displacement generator the quadruple (Ae_{1}, U, σ, V), where σ = (σ_{1}, …, σ_{k}). From the uniqueness of the SVD, the odr of a matrix is unique.

Observe that, since ||Δ(A)|| ≤ 2||Z|| · ||A|| ≤ 2||A|| and ||Δ(A)|| = σ_{1}, the orthogonal displacement representation of a matrix reduces the possibilities of cancellation in the computation of the matrix-vector product Ax.

We recall that multiplying an n × n matrix, represented through its odr (2.2) of length k, and a vector can be performed by means of 4k + 3 FFTs and 2k + 1 convolutions of length 2n.

The odr of a matrix A is strictly related to the ϵ-displacement representation. The following properties provide a characterization of the ϵ-displacement rank in terms of the SVD [4, 5].

Theorem 2.1. Let Δ(A) = UΣV^{T} be the singular values decomposition of Δ(A), where U, V ∈ ℝ^{n×k}, U^{T} U = I, V^{T} V = I, Σ = Diag(σ_{1}, …, σ_{k}) and σ_{1} ≥ … ≥ σ_{k} > 0 are the singular values of Δ(A). Let ϵ > be such that ϵ < σ_{1}. Then drk_{ϵ}(A) = r if and only if σ_{r} > ϵ ≥ σ_{r+1}. Moreover rdrk_{ϵ}(A) = r if and only if σ_{r} > ϵ||A|| ≥ σ_{r+1}.

Proof. We prove that σ_{r+1} ≤ ϵ if and only if drk_{ϵ}(A) ≤ r. This is sufficient to prove the theorem. Let u_{i}, v_{i}, i = 1, …, k, be the i-th column of U and V, respectively, that is

Δ(A) = ∑_{i=1}^{r} u_{i}σ_{i}v_{i}^{T} + ∑_{i=r+1}^{k} u_{i}σ_{i}v_{i}^{T} . (2.3)

Define E = ∑_{i=r+1}^{k} u_{i}σ_{i}v_{i}^{T}, then ||E|| = σ_{r+1}. Thus, if σ_{r+1} ≤ ϵ, we have drk_{ϵ}(A) ≤ drk(Δ(A)−E) = r. Vice-versa, if drk_{ϵ}(A) = r then there exists a matrix E, with ||E|| ≤ ϵ, such that rank(Δ(A)+E) = r. Define θ_{r} = min_{{X: rank(X)=r}} ||Δ(A)− X||; since θ_{r} = σ_{r+1} and σ_{r+1} = ||Δ(A) − ∑_{i=1}^{r} u_{i}σ_{i}v_{i}^{T} ||, it follows that σ_{r+1} = θ_{r} ≤ ||E||.

In particular, given 0 < ϵ < σ_{1} and the odg (Ae_{1}, U, σ, V), of A, it is possible to compute the ϵ-displacement representation of A of lower rank in the form

A_{ϵ} = L(Ae_{1}) + ∑_{i=1}^{r} σ_{i}L(u_{i})L(Zv_{i})^{T}

where r is such that σ_{r} > ϵ ≥ σ_{r+1}, or, in the case of the relative ϵ-displacement rank, σ_{r} > ϵ||A|| ≥ σ_{r+1}.

We call the above expression approximate orthogonal displacement representation (aodr) of A and the associated generator is called approximate orthogonal displacement generator (aodg).

The concept of ϵ-displacement rank has been implicitly used in [23] for implementing the Newton iteration for inverting Toeplitz matrices.

The following result provides a bound to the error in the ϵ-displacement representation of A.

Theorem 2.2. Let r = drk_{ϵ}(A) ≤ drk(A) = k, let A_{ϵ} be an ϵ-displacement representation of A, and let σ_{1}, …, σ_{k} be the singular values of δ(A). Then it holds

||A − A_{ϵ}|| ≤ n ∑_{i=r+1}^{k} σ_{i} ≤ n(k − r)ϵ.

Similarly, if r = rdrk_{ϵ}(A) and if A_{ϵ} is a relative ϵ-displacement representation of A, then it holds

\frac{||A_{ϵ} − A||}{||A||} ≤ n ∑_{i=r+1}^{k} σ_{i} ≤ n(k − r)ϵ.

Proof. Let Δ(A) = UΣV^{T} be the singular values decomposition of Δ(A). From Theorem 2.1 and (2.3) it holds A − A_{ϵ} = ∑_{i=r+1}^{k} σ_{i}L(u_{i})L(Zv_{i})^{T}, where ||u_{i}|| = ||v_{i}|| = 1. Since for any vector u it holds ||L(u)|| = ||L(u)^{T} || and ||L(u)|| ≤ ∑_{i=1}^{n} |u_{i}| · ||Z^{i}|| = ∑_{i=1}^{n} |u_{i}| ≤ \sqrt{n}||u||, we have ||A − A_{ϵ}|| ≤ n ∑_{i=r+1}^{k} σ_{i}. The proof is completed since |σ_{i}| ≤ ϵ for r = drk_{ϵ}(A) and |σ_{i}| ≤ ϵ||A|| for r = rdrk_{ϵ}(A), i = r + 1, …, k.

Conversely, we have the following

Theorem 2.3. If A has displacement rank k and B is an approximation of A such that ||A − B|| ≤ ϵ, then drk_{2ϵ}(B) ≤ k. Similarly, if ||B − A|| ≤ ϵ||A|| then rdrk_{2ϵ}(B) ≤ k.

Proof. It holds Δ(B) = Δ(A) + Δ(B − A), moreover, ||Δ(B − A)|| ≤ ||Z(B − A)|| + ||(B − A)Z|| ≤ 2||B − A|| ≤ 2ϵ since ||Z|| = 1.

The following inequalities are an immediate consequence of the definition:

drk_{δ}(A + B) ≤ drk_{ϵ}(A) + drk_{ϵ}(B),  δ = 2ϵ
drk_{δ}(AB) ≤ drk_{ϵ}(A) + drk_{ϵ}(B),  δ = ϵ(||A|| + ||B||)
drk_{δ}(A^{−1}) ≤ drk_{ϵ}(A),  δ = ||A^{−1}||ϵ
rdrk_{δ}(A + B) ≤ rdrk_{ϵ}(A) + rdrk_{ϵ}(B), δ = ϵ\frac{||A||+||B||}{||A+B||}
rdrk_{δ}(AB) ≤ rdrk_{ϵ}(A) + rdrk_{ϵ}(B),  δ = ϵ\frac{||A||·||B||}{||AB||}
rdrk_{δ}(A^{−1}) ≤ rdrk_{ϵ}(A),  δ = ϵμ(A),

(2.4)

where μ(A) = ||A|| · ||A^{−1}|| is the spectral condition number of ||A||.

Let us introduce the operator trunc_{ϵ}(·) defined on the sets of orthogonal displacement generators in the following way:

\begin{matrix}trunc_{ϵ}((a, U, σ, V)) = (a, Û, σ̂, V̂), \\ σ̂ = (σ_{1}, …, σk̂), Û = [u_{1}, …, uk̂], V̂ = [v_{1}, …, vk̂] ∈ ℝ^{n×k̂} \\ σk̂ > σ_{1}ϵ ≥ σ_{k̂+1}, \end{matrix} (2.5)

and with an abuse of notation, extend trunc_{ϵ}(·) to the set of matrices as A_{ϵ} = trunc_{ϵ}(A) where,

A_{ϵ} = L(Ae_{1}) + k̂ ; ∑ ; i=1 σ_{i}L(u_{i})L(ZJv_{i})^{T} .

Clearly from Theorem 2.2 it holds

||A − trunc_{ϵ}(A)|| ≤ n(k − k̂)σ_{1}ϵ

Equation (2.5) provides a means for computing an ϵ-displacement representation of A whenever an odg is available.

Remark 2.4. Constructing an orthogonal (ϵ)-displacement generator. Assume we are given F_{1}, F_{2} ∈ ℝ^{n×k} such that Δ(A) = F_{1}F_{2}^{T}. It is almost inexpensive to compute an orthogonal displacement generator of A, i.e. the SVD of Δ(A). A cheap way to compute the SVD of F_{1}F_{2}^{T} consists in computing the QR factorizations F_{1} = Q_{1}R_{1}, and F_{2} = Q_{2}R_{2}, so that we have Δ(A) = Q_{1}R_{1}R_{2}^{T} Q_{2}^{T}, and then the SVD R_{1}R_{2}^{T} = GΣH^{T} of the k × k matrix R_{1}R_{2}^{T}, where G^{T} G = I, H^{T} H = I.

Whence we obtain Δ(A) = UΣV^{T} where U = Q_{1}G, V = Q_{2}H. In this way, the overall cost is O(k^{2}n) ops. From the odg of A it is immediate to recover an ϵ-displacement generator by means of the operator trunc_{ϵ}(·).

2.1. A more convenient displacement operator. Multiplying a Toeplitzlike matrix and a vector by means of the odr of length k involves 4k + 3 FFTs and 2k + 1 convolutions of length 2n.

It is possible to reduce the lengths of FFTs and convolutions by a factor of 2 by means of a different displacement representation. The representation formulae that we introduce in this section are well known in the case of Toeplitz matrices [24, 23, 16, 13, 12, 8] and will be here derived, in their more general form, in terms of a suitable displacement operator.

Let us denote with C^{+} and C^{−} the matrices C^{+} = Z + e_{1}e_{n}^{T}, C^{−} = Z − e_{1}e_{n}^{T}, and set C^{+}(x) = ∑_{i=1}^{n} x_{i}(C^{+})^{i−1}, C^{−}(x) = ∑_{i=1}^{n} x_{i}(C^{−})^{i−1} the (+1)-circulant and (−1)-circulant matrices having x as first column. We recall that

C^{+}(x) = F Diag(y)F^{H},
C^{−}(x) = DF Diag(ŷ)F^{H}D,
y = \frac{1}{n}F^{H}x, ŷ = \frac{1}{n}F^{H}Dx,
F = (ω^{(i−1)(j−1)})_{i,j=1,... ,n}, ω = cos(2π/n) + i sin(2π/n)
D = Diag(1, θ, θ^{2}, …, θ^{n−1}), θ = cos(π/n) + i sin(π/n),

(2.6)

where i is the complex unit such that i^{2} = −1. Introduce the invertible operators

\begin{matrix}Δ^{+}(A) = C^{+}A − AC^{−} \\ Δ^{−}(A) = C^{−}A − AC^{+} \end{matrix}

such that

\begin{matrix}Δ^{+}(A) = −Δ(A) + e_{1}e_{n}^{T} A + Ae_{1}e_{n}^{T} \\ Δ^{−}(A) = −Δ(A) − e_{1}e_{n}^{T} A − Ae_{1}e_{n}^{T} . \end{matrix} (2.7)

We have the following result

Theorem 2.5. It holds

Δ^{+}(A) = ∑_{i=1}^{k} u_{i}v_{i}^{T} iff A = \frac{1}{2} ∑_{i=1}^{k} C^{+}(u_{i})C^{−}(Jv_{i}),
Δ^{−}(A) = ∑_{i=1}^{k} u_{i}v_{i}^{T} iff A = −\frac{1}{2} ∑_{i=1}^{k} C^{−}(u_{i})C^{+}(Jv_{i}),
Δ^{−}(A^{−1}) = −A^{−1}Δ^{+}(A)A^{−1},
Δ^{+}(AB) = Δ^{+}(A)B + AΔ^{+}(B) − 2Ae_{1}e_{n}^{T} B,
Δ^{−}(AB) = Δ^{−}(A)B + AΔ^{−}(B) + 2Ae_{1}e_{n}^{T} B,

where J is the permutation matrix having 1 on the anti-diagonal. In particular, if Δ^{+}(A) = ∑_{i=1}^{k} u_{i}v_{i}^{T} and det A ≠ 0, we obtain A^{−1} = \frac{1}{2} ∑_{i=1}^{k} C^{−}(A^{−1}u_{i})C^{+}(JA^{−T} v_{i}).

Proof. For the linearity of the operators Δ^{+} and Δ^{−} it is sufficient to prove the first two equations for k = 1, moreover, for the nonsingularity of the operators it is sufficient to prove the implications A = \frac{1}{2}C^{+}(u)C^{−}(Jv) ⇒ Δ^{+}(A) = uv^{T} and A = \frac{1}{2}C^{−}(u)C^{+}(Jv) ⇒ Δ^{−}(A) = uv^{T}. For A = \frac{1}{2}C^{+}(u)C^{−}(Jv) we have Δ^{+}(A) = C^{+}A − AC^{−} = \frac{1}{2}(C^{+}(u)C^{+}C^{−}(Jv) − C^{+}(u)C^{−}C^{−}(Jv)) = C^{+}(u)e_{1}e_{n}^{T} C^{−}(Jv) = uv^{T}. Similarly we may do for the operator Δ^{−}. The remaining equations are straightforward.

Observe that, for a Toeplitz matrix A it holds Δ^{+}(A) = e_{1}(u − v)^{T} J + (u + v)e_{1}^{T} J, for u = Ae_{1}, v = ZAe_{n}, that is, A = \frac{1}{2}(C^{+}(u + v) + C^{−}(u − v)), therefore we have the following inversion formula (compare with [24, 23, 16, 13, 12, 8]) A^{−1} = −\frac{1}{2}(C^{−}(A^{−1}e_{1}) (I − C^{+}(w)) + (I + C^{−}(w)) C^{+}(A^{−1}e_{1})), w = A^{−1}ZAe_{n}.

The new displacement representation given in the above theorem allows one to compute the product Ax by means of 4k + 2 FFT’s, 2k convolutions and 2k + 1 scalings of length n.

Indeed, we may apply the concepts of orthogonal displacement representation and of ϵ-displacement rank to the representations obtained through the new operators Δ^{+} and Δ^{−}. In this way we arrive at a more efficient implementation of this tool.

Due to the nonsingularity of Δ^{−} and Δ^{+}, the odg associated with Δ^{−} and Δ^{+}, is simply defined by the triple (U, σ, V) such that UΣV is the SVD of Δ^{−}(A) or of Δ^{+}(A), respectively, where Σ = Diag(σ_{1}, …, σ_{k}).

Observe also that the displacement ranks obtained by means of the operators Δ, Δ^{+} and Δ^{−} may differ at most by 2 and from the relations (2.1,2.2,2.7) we obtain the representation formula

\begin{matrix}A = \frac{1}{2}(∑_{i=1}^{k} C^{+}(u_{i})C^{−}(Jv_{i}) + C^{−}(JA^{T} e_{n}) + C^{+}(Ae_{1})) \\ Δ(A) = − ∑_{i=1}^{k} u_{i}v_{i}^{T} . \end{matrix} (2.8)

2.2. Matrices with additional structure. Further reduction of the computational costs can be obtained if the matrix A has more structure. For instance, if A is Toeplitz, then it is symmetric with respect to the secondary diagonal (persymmetric), and it holds A^{T} = JAJ, where J is the permutation matrix having 1 on the secondary diagonal. The persymmetry of A implies that Δ(A)J is skew-symmetric.

For skew-symmetric matrices we have the following technical result:

Lemma 2.6. Let A ∈ ℝ^{n×n} be such that A = −A^{T} and rank(A) = k. Then k = 2h is even and there exist W ∈ ℝ^{n×k}, D ∈ ℝ^{h×h} such that W^{T} W = I_{k}, D = Diag(d_{1}, …, d_{h}) and

A = W (D ⊗ S)W^{T} , S = [ \begin{matrix}0 & 1 \\ −1 & 0 \end{matrix} ] ,

moreover ±id_{j} are the nonnull eigenvalues of A, where i is the complex unit such that i^{2} = −1, and |d_{j}| = σ_{2j−1} = σ_{2j}, j = 1, …, h, where σ_{1}, …, σ_{k} are the singular values of A.

Proof. From the Schur normal form deduce that a skew symmetric matrix A is such that A = QDQ^{H}, for Q^{H}Q = I, D = Diag(d_{j}) where d_{2j−1} = θ_{j}i, d_{2j} = −θ_{j}i. Therefore, only an even number of eigenvalues can be nonzero. Denote q_{j} = Qe_{j}, where e_{j} denotes the j-th column of the identity matrix, apply the equation

[ \begin{matrix}θi & 0 \\ 0 & −θi \end{matrix} ] = \frac{θ}{2} [ \begin{matrix}1 & −i \\ 1 & i \end{matrix} ] [ \begin{matrix}0 & 1 \\ −1 & 0 \end{matrix} ] [ \begin{matrix}1 & 1 \\ i & −i \end{matrix} ]

and conclude that the claim is true where the columns of W are given by (q_{2j−1} + q̄_{2j−1})\sqrt{2} and (q_{2j−1} − q̄_{2j−1})i\sqrt{2}, that are orthonormal.

As a consequence, we have the following fundamental property:

Theorem 2.7. [Displacement representation of a persymmetric matrix]. If the persymmetric matrix A has displacement rank k with respect to Δ, then k = 2h and there exist matrices W ∈ ℝ^{n×k}, D = Diag(d_{1}, …, d_{h}) ∈ ℝ^{n×n} such that |d_{1}| ≥ |d_{2}| ≥ … ≥ |d_{h}|, W^{T} W = I_{k} and

Δ(A) = W (D ⊗ S)W^{T} J, S = [ \begin{matrix}0 & 1 \\ −1 & 0 \end{matrix} ] ,

Therefore,

A = L(a) + ∑_{i=1}^{h} d_{i}(L(w_{2i−1})L^{T} (ZJw_{2i}) − L(w_{2i})L^{T} (ZJw_{2i−1})) (2.9)

where a = Ae_{1}, w_{i} = W e_{i}.

Operating with the persymmetric odr (2.9) is more convenient in terms of memory space and computational cost. For instance, the number of FFT’s sufficient for computing the product of A and a vector can be reduced to 3k + 3.

For the operators Δ^{+} and Δ^{−} the persymmetry property provides a more structured representation that can be recovered from (2.8) where A^{T} e_{n} = JAe_{1}, in the light of Theorem 2.7. Also in this case the number of FFTs sufficient to compute the matrix-vector product is reduced to 3k + 3.

3. Applications

The displacement representation allows one to compute the product between a matrix with displacement rank k and a vector in O(kn log n) operations by means of FFT. Similarly the ϵ-displacement representation allows one to approximate the matrix-vector product with the same asymptotic cost.

In this section we present two applications where sequences of matrices are generated with a displacement rank that grows exponentially while their ϵ-rank remains bounded still guaranteeing the convergence properties of the original algorithm and at the same time, reducing the complexity of the computation.

The two applications concern the inversion of Toeplitz-like matrices by means of Newton’s iteration, and the solution of block banded block Toeplitz system with banded Toeplitz blocks (BBBT-BTB) with cyclic reduction.

3.1. Newton’s Iteration. Newton’s iteration is a valid tool for solving problems where the solution to be computed can be formally viewed as the zero of a known function. This method has been analyzed in [8] and, concerning Toeplitz matrices, has been investigated in [23] where the iteration (1.2) has been modified by replacing each X_{i} with the closest matrix having displacement rank 2. Condition on X_{0} have been given in order that the method, modified in this way, still converges.

Here we present a different approach, based on the concept of adaptive algorithm, where each approximation X_{i} is replaced by the closest matrix having displacement rank k_{i}, and k_{i} is chosen so that the distance of A from the matrices with displacement rank k is “sufficiently” small in order to guarantee the convergence.

In this way, the displacement rank is not bounded a priori but it may vary according to the difficulty of the problem and to the needs of the computation. This approach is adaptive in the sense that the cost of the algorithm is not uniform with respect to the set of input data, where the algorithm adjusts the computation to the specific difficulty of the problem. In practice, for well conditioned matrices the growth of the ϵ-displacement rank is negligible and the overall cost is O(n log n) ops.

For (polynomially) ill-conditioned matrices the cost seems to grow as O(n log^{2} n) ops. In fact, our modification of Newton’s iteration requires O(kn log n) ops per step, where k is the maximum displacement rank of the matrices X_{i} involved in the computation. In the numerical experiments that we report at the end of this section, k takes very moderate values (say, in the range [2,6]) even for very ill-conditioned matrices.

Let A be a nonsingular matrix and consider the sequence of matrices generated, from an assigned X_{0}, by means of (1.2). This sequence that converges to A^{−1} if det A ≠ 0 for suitable X_{0}, is formally obtained by applying Newton’s iteration x_{i+1} = x_{i} − f(x_{i})/f′(x_{i}) to the function f(X) = X^{−1} − A. Moreover we may look at (1.2) as to a particular instance of the generalized Newton’s iteration

X_{i+1} = A^{−1}(I − (I − AX_{i})^{m}), m ≥ 2, (3.1)

which, for m = 2 reduces to (1.2) and for m = 3 to

X_{i+1} = 3X_{i} − 3X_{i}AX_{i} + X_{i}(AX_{i})^{2}. (3.2)

3.1.1. Convergence and regularization properties. Denoting R_{i} = I − AX_{i} the residual error at the i-th step, from (3.1) we obtain that R_{i+1} = R_{i}^{m}, that is,

R_{i} = R_{0}^{m^{i}}.

The above relation provides a lot of information on the convergence of the general iteration (3.1).

For positive definite A, with eigenvalues 0 < λ_{1} ≤ λ_{2} ≤ … ≤ λ_{n} = ρ(A), and with X_{0} = (1/θ)I, θ ≥ ρ(A), we find that R_{0} has eigenvalues 1 − λ_{i}/θ between 0 and 1. Therefore ||R_{i}|| = (1 − λ_{1}/θ)^{m^{i}} gives the double exponential convergence since 0 ≤ 1 − λ_{1}/θ ≤ 1 − 1/μ(A) < 1. A simple choice for θ can be θ = ||A|| for any induced norm || · ||.

For a general A the choice X_{0} = (1/θ)A^{T}, where θ ≥ ρ(A^{T} A), guarantees the double exponential convergence since R_{0} = I − (1/θ)AA^{T} is positive definite and has eigenvalues less than or equal to 1 − μ_{1}/θ, where μ_{1} is the smallest eigenvalue of AA^{T}.

Observe that, since (1 − 1/x)^{x} < 1/e < 1/2, after k = ⌈log_{m}(λ_{1}/θ)⌉ iterations, we have ρ(R_{k}) ≤ 1/e < 1/2, therefore, in order to have the residual error less than a given error bound ϵ, it is sufficient to apply the iteration (3.1) for ν = ⌈log_{m}(λ_{1}/θ)⌉ + log_{m} log ϵ^{−1} steps, where μ(A) ≤ λ_{1}/θ.

The first additive term in the above expression, that we call overhead number of iterations, provides the number of steps needed in order to arrive at a pure double exponential convergence of order m, where for pure double exponential convergence of order m we mean the error tends to zero as a^{m^{i}} for a ≤ 1/2 and m > 1.

The second term provides the number of steps of pure double exponential convergence that occurs to reach the residual error bound ϵ.

The number of overhead iterations grows roughly as log μ(A) and for Toeplitz matrices associate with a symbol f(x) > 0 the condition number is independent of n, while for f(x) ≥ 0 the condition number grows at most polynomially with n, therefore the overhead number of iteration is O(log n).

Better choices for X_{0} can be done by setting X_{0} = (1/θ)P, if A is positive definite, or X_{0} = (1/θ)P P^{T} A^{T} where P is a preconditioner of A, i.e., a matrix such that μ(AP), is smaller than μ(A) and P has small displacement rank. There is a wide literature concerning such preconditioners for the conjugate gradient method, that can be exported to the case of Newton’s iteration. For instance, the least square circulant preconditioner of [11] provides an effective reduction of the condition number. Moreover, the computation of P is quite easy, since P = C^{−1} where the first column c = (c_{0}, …, c_{n−1}) of the circulant matrix C is given by c_{0} = a_{0}, c_{i} = ((n − i)a_{i} + ia_{n−i})/n, i = 1, …, n − 1, where a_{0}, a_{1}, … a_{n−1} are the entries in the first column of A.

If P is chosen as the inverse of the banded Toeplitz matrix whose associated function is a trigonometric polynomial interpolating f(z) and its derivatives at its roots, then μ(CA) is bounded from above by a constant [9, 10].

Newton’s iteration has also very nice regularizing properties.

Remark 3.1. On the regularizing properties of Newton’s iteration. Consider for simplicity the case where A is positive definite. Denoting q_{j}, j = 1, …, n the eigenvectors of A corresponding to the eigenvalues λ_{1}, …, λ_{n}, from the relation A = ∑_{j=1}^{n} λ_{j}q_{j}q_{j}^{T}, we deduce that R_{i} = ∑_{j=1}^{n}(1− \frac{λ_{j}}{θ})^{m^{i}}q_{j}q_{j}^{T}. That is, the components of the residual error along the matrices q_{j}q_{j}^{T} converge to zero much faster for j close to n than for j close to 1, that is, for the eigenvectors q_{j} corresponding to the largest eigenvalues. Similarly we have

A^{−1} − X_{i} = ∑_{j=1}^{n} \frac{(1 − λ_{j}/θ)^{m^{i}}}{λ_{j}} q_{j}q_{j}^{T} .

For very ill-conditioned matrices, X_{i} provides a poor approximation to A^{−1} in the eigenspace corresponding to the smallest eigenvalues and a very good approximation in the eigenspace corresponding to the largest eigenvalues even for moderate values of i. This property is very important since it provides a sort of “regularizing” behavior of Newton’s iteration in the sense that a few iteration of (1.2) provide a generalized inverse of A that does not amplify the possible error in the known term of the linear system since the smallest eigenvalues of A are not yet inverted. This feature, typical of the conjugate gradient technique is called regularization.

3.1.2. Using ϵ-displacement rank: adaptive Newton’s iteration. Indeed the Newton iteration for general matrices is very expensive since each step has a cost O(n^{3}) ops versus O(n^{2}) ops of the PCG. However, Newton’s iteration may be competitive in the case of certain structured matrices where the same structure of A is shared by A^{−1}. This is the case of the displacement structure.

If for a Toeplitz (or more generally a Toeplitz-like) matrix we are able to implement a single Newton’s iteration in O(n log n) steps, then we obtain an adaptive solver that require a number of ops ranging from O(n log n log log ϵ^{−1}) to O(n log^{2} n + n log n log log ϵ^{−1}) according to the conditioning of the matrix A.

Let drk(A) = k; we will analyze the modification of Newton’s iteration defined by the following recurrence

\begin{matrix}Y_{i} = 2X_{i} − X_{i}AX_{i} \\ X_{i+1} = trunc_{ϵ_{i}}(Y_{i}) \end{matrix} (3.3)

where trunc_{ϵ}(·) is defined in (2.5) and the choice of ϵ_{i} can be made according to different strategies. The simplest one is to choose a constant value of ϵ_{i}, say, the square root of the machine precision of the computation.

A dynamic choice of ϵ that guarantees the convergence of the modified Newton method can be easily designed. Before delving in this problem observe that

Δ(Y_{i}) = 2Δ(X_{i}) − Δ(X_{i})AX_{i} − X_{i}Δ(A)X_{i} − X_{i}AΔ(X_{i})

Δ^{−}(Y_{i}) = 2Δ^{−}(X_{i}) − Δ^{−}(X_{i})AX_{i} − X_{i}Δ^{+}(A)X_{i} − X_{i}AΔ^{−}(X_{i}) (3.4)

so that drk(Y_{i}) ≤ k + drk(X_{i}) both for Δ and Δ^{−}, where k = drk(A).

Moreover, replacing for notational simplicity X_{i} with X and Y_{i} with Y, if Δ(X) = U_{X}D_{X}V_{X}^{T}, or Δ^{−}(X) = U_{X}D_{X}V_{X}^{T}, where U_{X}, V_{X} ∈ ℝ^{n×h} are orthogonal, and Δ(A) = U_{A}D_{A}V_{A}^{T}, or Δ^{+}(A) = U_{A}D_{A}V_{A}^{T}, where U_{A}, V_{A} ∈ R^{n×k} are orthogonal, we obtain that

\begin{matrix}Δ(Y) = U_{Y} D_{Y} V_{Y}^{T} \\ U_{Y} = [U_{X}, XU_{A}, XAU_{X}], V_{Y} = [V_{X}, X^{T} U_{A}, X^{T} A^{T} U_{X}] \\ D_{Y} = [ \begin{matrix}2D_{X} & −D_{X} & 0 \\ −D_{X} & 0 & 0 \\ 0 & 0 & −D_{A} \end{matrix} ] \end{matrix} (3.5)

Thus, an odr of X_{i+1} = trunc_{ϵ}(Y) can be computed by performing the following steps:

Algorithm 1. Computing the aodr of Y = 2X − XAX (operator Δ^{−}).

1. Compute the matrices U_{Y}, V_{Y}, D_{Y} of (3.5);

2. Compute the QR factorizations U_{Y} = Q_{1}R_{1}, V_{Y} = Q_{2}R_{2};

3. Compute the SVD R_{1}D_{Y} R_{2}^{T} = UΣV, set Σ̂ = Diag(σ_{1}, …, σ_{r}), where r is such that for the singular values of R_{1}D_{Y} R_{2}^{T} it holds σ_{r+1} ≤ ϵσ_{1} < σ_{r};

4. Set U_{X} = Q_{1}Û, V_{X} = Q_{2}V̂, D_{X} = Σ̂, where Û and V̂ are made up by the first r columns of U and V, respectively.

The above algorithms can be extended to the operator Δ; in this case we need also to compute the first column of X_{i+1}, i.e., 2X_{i}e_{1}−X_{i}AX_{i}e_{1} in order to complete the odg.

The cost of the above computation is about O(hk) FFTs. Improvements of this computation can be obtained by suitably using further structures, like the persymmetry (as in the case of Toeplitz matrices) or the symmetry of the matrices involved, that allow the reduction of the number of FFTs. These issues, more closely related to the efficient implementation of the algorithm, are not treated in this paper.

Concerning the design of strategies for the determination of the cutting level ϵ, let us denote D_{X} = (σ_{1}, …, σ_{h}), Y = 2X − XAX, V = trunc_{ϵ}(Y) = Y + E, R_{X} = I − AX, R_{Y} = I − AY, R_{V} = I − AV.

It clearly holds drk(Y) ≤ 2h + k, R_{Y} = R_{X}^{m}, ||E|| ≤ 2n ∑_{i=r+1}^{2h+k} σ_{i} ≤ 2n(2h + k − r)ϵσ_{1} where σ_{r} > ϵσ_{1} ≥ σ_{r+1}, and R_{V} = R_{Y} − AE. Therefore ||R_{V} || ≤ ||R_{X}^{m}|| + ||A||2n(2h + k − r)ϵ. Thus, we obtain that

||R_{V} || ≤ ||R_{X}^{m}||(1 + δ), for ϵ ≤ δ||R_{X}^{m}||/(||A||2n(2h + k)σ_{1}). (3.6)

The above inequality is the basis to prove the following

Theorem 3.2. Let X_{0} be such that ||X_{0}|| ≤ 1 − θ, set ϵ_{i} = θ/(4n(2h_{i} + k)σ_{1}^{(i)}||A||), h_{i} = drk(X_{i}), σ_{1}^{(i)} = ρ(Δ(X_{i})) (σ_{1}^{(i)} = ρ(Δ^{−}(X_{i}))). Then for the residuals R_{i} = I − AX_{i} of the sequence {X_{i}}_{i} generated by the modified Newton iteration it holds ||R_{i}|| ≤ (1 − θ/2)^{m^{i}}.

Proof. From 3.6 it holds ||R_{i}|| ≤ ||R_{i−1}||^{m}(1 + θ/2). Whence we obtain ||R_{i}|| ≤ (1 − θ)(1 + θ/2)^{m^{i}} ≤ (1 − θ/2)^{m^{i}}.

For θ = 1/μ(A), the bound on ϵ turns into ϵ ≤ ||R_{i}||^{m}/(4n(2h_{i} + k)σ_{1}μ(A)). Indeed, this bound guarantees the convergence whatever is the matrix A. However, in practice is too pessimistic and may lead to an unnecessary large growth of the displacement rank of the intermediate X_{i}. A more realistic strategy is to choose ϵ according to specific heuristics like, for instance, ϵ_{i} = max\{||R_{i}||/μ(A), \sqrt{u}\}, where u is the machine precision of the floating point arithmetic used. The tables 1,2, 3, 4 report the values of the displacement and of the residual at the intermediate steps of the modified Newton’s iteration, with no preconditioning, for different classes of matrices.

We observe that the growth of the displacement rank depends on the condition number of the matrix. However, it is much inferior to the size n of the matrix in all the test performed.

Below we resume our modification of Newton’s iteration for the inversion of a Toeplitz-like matrix assigned in terms of its odr with respect to the operator Δ^{−}. A similar formulation can be given for the operator Δ.

Algorithm 2. Newton’s iteration with ϵ-displacement rank (operator Δ^{−})

Input: Integers n, k, the maximum number MAXIT of iterations, a residual error bound ϵ, the odg (U_{A}, σ_{A}, V_{A}) of the nonsingular matrix A ∈ ℝ^{n×n} having displacement rank k; the odg (U_{X}, σ_{X}, V_{X}) of a matrix X_{0} such that AX_{0} is positive definite and ρ(AX_{0}) < 1; a strategy for the choice of the truncation value ϵ_{i}.

Output: An approximation X of A^{−1} such that ||I − AX|| ≤ ϵ given in terms

of its odg (U_{X}, σ_{X}, V_{X}). Computation:

Stage 0: Set i = 0, X = X_{0};

Stage 1: Compute the truncation value ϵ_{i};

Stage 2: Compute the aodg (U_{Y} , σ_{Y} , V_{Y}) of Y = 2X−XAX, by means of algorithm 1 with ϵ = ϵ_{i} and set (U_{X}, σ_{X}, V_{X}) = (U_{Y} , σ_{Y} , V_{Y}) Stage 3: Compute ||R|| for the residual R = I − AX; if ||R|| > ϵ and

i ≤ MAXIT set i = i + 1 and continue from Stage 1. Stage 4: If i > MAXIT then output FAILURE, otherwise output

odg(X).

Concerning the choice of ϵ_{i}, Theorem 3.2 provides the strategy that guarantees convergence. Different strategies based on heuristics can be applied. The starting matrix X_{0} can be chosen according to the features of A: for a general A, we may set X_{0} = A^{T} /||AA^{T} ||; for A symmetric and positive definite we may set X_{0} = I/||A||. A more convenient choice can be done in terms of the least square circulant approximating A.

The computation of ||V ||, for a matrix V assigned in terms of its odg of length h, is generally expensive since it requires the computation of all the n^{2} entries of V. An alternative way is to replace ||V || with an upper bound obtained from its odr, say, ||V ||_{1} ≤ n ∑_{i=1}^{h} σ_{i}. This bound can be used for estimating ||A|| as well as ||A^{T} A||, and for computing ||R|| at stage 3. The latter estimate requires the

n Nstep Mdrk Sdrk err
32  8  3  16 2.4e-14
64  8  3  16 5.1e-14
128  8  3  16 6.6e-14
256  8  3  16 6.8e-14
512  8  3  16 6.0e-14
1024  8  3  16 5.8e-14
2048  8  3  16 5.3e-14
4096  8  3  16 4.7e-14

Table 1. Matrix 1

n Nstep Mdrk Sdrk err
32  8  3  15 2.0e-16
64  8  3  15 4.0e-16
128  8  3  15 2.0e-16
256  8  3  15 2.0e-16
512  8  3  15 2.0e-16
1024  8  3  15 2.1e-16
2048  8  3  15 2.0e-16
4096  8  3  15 2.2e-16

Table 2. Matrix 2

computation of the odg of R. A cheaper way that provides a necessary but, in principle, not sufficient condition in order that ||R|| ≤ ϵ is ||Re_{1}|| ≤ ϵ.

3.1.3. Numerical experiments. In Tables 1, 2, 3, 4 we report the number of steps Nsteps, the maximum displacement rank Mdrk, the sum of the displacement ranks Sdrk of each step, and the residual error err for inverting the matrix A_{n}, for each value of n, where A_{n} = (a_{|i−j|})_{i,j=1,n} is defined in the following way:

Matrix 1: a_{|i−j|} = 1/(|i − j| + 1);
Matrix 2: a_{0} = 4, a_{1} = 1, a_{i} = 0 otherwise;
Matrix 3: a_{0} = 2.1, a_{1} = 1, a_{i} = 0 otherwise;
Matrix 4: a_{0} = 2, a_{1} = 1, a_{i} = 0 otherwise;

The quantity Sdrk · n log n expesses the overall arithmetic cost of the algorithm.

The matrices 1,2,3 are asymptotically well conditioned while the matrix 4 has condition number that grows as O(n^{2}).

For well conditioned matrices the number of iterations as well as the maximum displacement rank is independent of n. For the matrix 4, the number of iterations as well as the sum of the displacement ranks seems to grow as log_{2} n.

The number of iterations, that grows proportionally to log μ(A), can be reduced to a constant by using a band matrix as preconditioner.

n Nstep Mdrk Sdrk err
32  11  4  38 3.3e-14
64  12  5  41 1.7e-15
128  11  8  52 2.0e-15
256  11  7  47 2.1e-15
512  10  6  41 3.4e-13
1024 10  5  38 2.5e-13
2048 10  5  35 1.9e-13
4096 10  4  29 2.2e-13

Table 3. Matrix 3

n Nstep Mdrk Sdrk err
32  12  5  47 3.8e-15
64  13  7  57 2.5e-15
128  14  8  68 8.9e-15
256  15  8  76 4.8e-14
512  16  10  93 6.5e-14
1024 18  13  121 1.1e-13
2048 22  19  187 3.2e-13
4096 23  22  224 1.3e-13

Table 4. Matrix 4

3.2. Cyclic reduction. The solution of block-tridiagonal block-Toeplitz systems

[ \begin{matrix}A & B \\ C & A & ... \\ & ... & ... & B \\ & & C & A \end{matrix} ] [ x_{1} ; x_{2} ; ⋮ ; x_{n} ] = [ b_{1} ; b_{2} ; ⋮ ; b_{n} ] (3.7)

of block-size n = 2^{q} +1, can be efficiently performed by means of cyclic reduction [1, 2, 3, 25]. In this computation the sequence of matrices A_{i}, B_{i}, C_{i}, i = 1, 2, …, q, such that

\begin{matrix}A_{i+1} = A_{i} − B_{i}A_{i}^{−1}C_{i} − C_{i}A_{i}^{−1}B_{i} \\ B_{i+1} = −B_{i}A_{i}^{−1}B_{i} \\ C_{i+1} = −C_{i}A_{i}^{−1}C_{i} \\ A_{0} = A, B_{0} = B, C_{0} = C, \end{matrix} (3.8)

must be computed, and the block components x_{j}, j = 1, …, n can be easily recovered from the blocks A_{i}, B_{i}, C_{i} and the vectors b_{j}.

In certain computations, for instance when the system (3.7) is obtained by reblocking a banded Toeplitz matrix [3] or in queuing theory problems [6], the blocks (3.8) have displacement rank bounded by a small constant.

When the system (3.7) is obtained by reblocking a block-banded block-Toeplitz system with banded Toeplitz blocks, the matrices (3.8) have displacement rank that may grow exponentially until it reaches the size of the blocks. This is typical in image processing problems. However, in [4, 5] it is observed that the approximate displacement rank of the matrices (3.8) is bounded by a small constant. Therefore a reformulation of the cyclic reduction in terms of aodr can be given.

We refer the reader to [4] for details on cyclic reduction where the blocks (3.8) are approximated by aodr. Here we describe how the aodg of the blocks (3.8) can be computed.

Observe that from the properties of the operators Δ and Δ^{+} we have

Δ(BA^{−1}B) = Δ(B)A^{−1}B − BA^{−1}Δ(A)A^{−1}B + BA^{−1}Δ(B)

Δ^{+}(BA^{−1}B) = Δ^{+}(B)A^{−1}B − BA^{−1}Δ^{+}(A)A^{−1}B + BA^{−1}Δ^{+}(B).

Similar relations hold for CA^{−1}C and A − BA^{−1}C − CA^{−1}B.

The above formulae can be used in the same way as (3.4,3.5) for the design of an algorithm, like Algorithm 1 for the computation of one step of approximate cyclic reduction. The computational cost of the algorithm that we obtain for the approximate solution of an n × n block-banded block-Toeplitz system with banded Toeplitz blocks, where n = m^{2} and the blocks have size m blocks, is O(h^{2}m log^{2} m log n + hnm log m) arithmetic operations, where h is the largest ϵdisplacement rank encountered in the overall computation.

The algorithm has been implemented and tested for a problem of image restoration of N × N pixels images, with the PSF defined by the mask

σ [ \begin{matrix}1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & α & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 \end{matrix} ] , σ = 1/(24 + α).

This mask generates a 5 block-diagonal block-Toeplitz matrix of block-size N, with blocks of size N; the off diagonal blocks are 5-diagonal Toeplitz matrices with all the nonzero entries equal to σ, the diagonal blocks are 5-diagonal Toeplitz matrices having σα on the main diagonal and σ on the remaining nonzero diagonals. The matrix obtained in this way is reblocked into a block tridiagonal form, thus obtaining (3.7).

We have considered images of size N = 64, that generate a block tridiagonal matrix with block size n = 32, with blocks of size m = 128. We have reconstructed the image by solving the system with our modified CR algorithm. For α = 5 the conditions that guarantee the boundness of the approximate displacement rank are satisfied. We have reconstructed the image with values of ϵ in the range [10^{−6}, 10^{−3}], keeping a small value of the ϵ-displacement rank and an accurate precision in the approximation. Larger values of ϵ did not provide adequate approximation of the solution. The original, blurred and restored images are reported in figures 1, 2.

We have also computed the maximum values of drk_{ϵ} and rdrk_{ϵ} of the blocks generated by CR for different values of α, ϵ and the original block size N. The results are summarized in tables 5, 6, where we report the size m of the blocks of the reblocked tridiagonal matrix and the largest values of drk_{ϵ} and rdrk_{ϵ} of the blocks generated by CR, for different values of ϵ. We observe that for α = 2 the hypotheses that guarantee the conservation of the ϵ-displacement rank are not satisfied [4] and even though the rdrk_{ϵ} depends on N, its growth is rather moderate.

16 DARIO ANDREA BINI AND BEATRICE MEINI

Figure 1. Original and blurred images Figure 2. Restored images, with ϵ = 10^{−3} and ϵ = 10^{−6}

ϵ = 10^{−2} ϵ = 10^{−4} ϵ = 10^{−6}
m drk_{ϵ} rdrk_{ϵ} drk_{ϵ} rdrk_{ϵ} drk_{ϵ} rdrk_{ϵ}
64 12  12  20  22  28  34
128 12  12  20  22  28  34
256 12  12  20  22  28  34

Table 5. α = 8

\begin{matrix}m & drk_{ϵ} & rdrk_{ϵ} \\ 64 & 52 & 22 \\ 128 & 88 & 42 \\ 256 & 150 & 62 \end{matrix}

Table 6. α = 2, ϵ = 10^{−2}

References

1. D. A. Bini and B. Meini, On cyclic reduction applied to a class of Toeplitz-like matrices arising in queueing problems, Computations with Markov Chains (W. J. Stewart, ed.), Kluwer Academic Publisher, 1995, pp. 21–38.

2. , Improved cyclic reduction for solving queueing problems, Numerical Algorithms 15 (1997), 57–74.

3. , Effective methods for solving banded Toeplitz systems, SIAM J. Matrix Anal. Appl. 20 (1999), 700–719.

4. , Solving block banded block Toeplitz systems with banded Toeplitz blocks, Advanced Signal Processing Algorithms, Architectures, and Implementations IX (F. T. Luk, ed.), SPIE-The International Society for Optical Engineering, 1999, In press.

5. , Solving block banded block Toeplitz systems with structured blocks: new algorithms and open problems, Proceedings of WLSSC 1999, Sozopol, Bulgaria (P.Y. Yalamov S. D. Margenov, P. S. Vassilevski, ed.), 1999, In press.

6. D. A. Bini, B. Meini, and V. Ramaswami, Analyzing M/G/1 paradigms through QBDs: the role of the block structure in computing the matrix G, submitted for publication, 1999.

7. D. Bini and V. Pan, Improved parallel computations with Toeplitz-like and Hankel matrices, Linear Algebra Appl. 188–189 (1993), 3–29.

8. D. A. Bini and V. Pan, Matrix and polynomial computations, vol. 1: Fundamental algorithms, Birkhäuser, Boston, 1994.

9. R. H. Chan, Toeplitz preconditioners for Toeplitz systems with nonnegative generating functions, IMA J. of Numerical Analysis 11 (1991), 333–345.

10. R. H. Chan and M. K. Ng, Iterative methods for linear systems with matrix structure, Fast reliable algorithms for matrices with structure (T. Kailath and A. H. Sayed, eds.), SIAM, Philadelphia, 1999, pp. 117–154.

11. T. Chan, An optimal circulant preconditioner for Toeplitz systems, SIAM J. Sci. Stat. Comput. 9 (1988), 766–771.

12. I. C. Gohberg, T. Kailath, and V. Olshevsky, Fast Gaussian elimination with partial pivoting for matrices with displacement structure, Math. Comp. 64 (1995), 1557–1576.

13. G. Heinig, Inversion of generalized Cauchy matrices and other classes of structured matrices, Linear Algebra, Signal Process., IMA Vol. Math. Appl. 69 (1995), 95–114.

14. , Transformation approaches for fast and stable solutions of Toeplitz systems and polynomial equations, Proceedings of the international workshop ”Recent advances in applied mathemetics”, State of Kuwait, May 4–7, 1996, pp. 223–238.

15. G. Heinig and K. Rost, Algebraic methods for Toeplitz-like matrices and operators, Akademie-Verlag, Berlin, and Birkhäuser, Boston, 1984.

16. T. Kailath, Displacement structure and array algorithms, Fast reliable algorithms for matrices with structure (T. Kailath and A. H. Sayed, eds.), SIAM, Philadelphia, 1999, pp. 1–56.

17. T. Kailath, S. Kung, and M. Morf, Displacement rank of a matrix, Bulletin of the American Mathematical Society 1 (1979), 769–773.

18. , Displacement ranks of matrices and linear equations, J. Math. Anal. Appl. 68 (1979), 395–407.

19. T. Kailath and V. Olshevsky, Displacement structure approach to Chebyshev-Vandermonde and related matrices, Integral Equations Operator Theory 22 (1995), 65–92.

20. , Displacement structure approach to discete-trigonometric-transform based preconditioners of G. Strang type and T. Chan type, Calcolo 33 (1996), 191–208.

21. T. Kailath and A. H. Sayed, Displacement structure: theory and applications, SIAM Review 37 (1995), 297–386.

22. T. Kailath and A. H. Sayed (eds.), Fast reliable algorithms for matrices with structure, SIAM, Philadelphia, 1999.

23. V. Y. Pan, S. Branham, R. E. Rosholt, and A.-L. Zheng, Newton’s iteration for structured matrices, Fast reliable algorithms for matrices with structure (T. Kailath and A. H. Sayed, eds.), SIAM, Philadelphia, 1999, pp. 189–210.

24. M. Van Barel, G. Heinig, and P. Kravanja, A stabilized superfast solver for nonsymmetric Toeplitz systems, Tech. Report TW 293, Katholieke Universiteit Leuven, Department of Computer Science, October 1999.

25. P. Yalamov and V. Pavlov, Stability of the block cyclic reduction, Linear Algebra Appl. 249 (1996), 341–358.

Dipartimento di Matematica, Università di Pisa, via Buonarroti 2, I-56127 Pisa, Italy

E-mail address: bini@dm.unipi.it, meini@dm.unipi.it