On the regularized solution of block banded block Toeplitz systems

D. A. Bini^{a}, A. Farusi^{b}, G. Fiorentino^{a} and B. Meini^{a}

aDipartimento di Matematica, via Buonarroti 2, 56127 Pisa, Italy bIstituto di Elaborazione dell’Informazione, CNR, S. Cataldo, Pisa, Italy

ABSTRACT

Given two n × n Toeplitz matrices T_{1} and T_{2}, and a vector b ∈ R^{n^{2}}, consider the linear system Ax = b − η, where η ∈ R^{n^{2}} is an unknown vector representing the noise and A = T_{1} ⊗ T_{2}. Recovering approximations of x, given A and b, is encountered in image restoration problems. We propose a method for the approximation of the solution x that has good regularization properties. The algorithm is based on a modified version of Newton’s iteration for matrix inversion and relies on the concept of approximate displacement rank. We provide a formal description of the regularization properties of Newton’s iteration in terms of filters and determine bounds to the number of iterations that guarantee regularization. The method is extended to deal with more general systems where A = ∑_{i=1}^{h} T_{1}^{(i)}⊗T_{2}^{(i)}. The cost of computing regularized inverses is O(n log n) operations (ops), the cost of solving the system Ax = b is O(n^{2} log n) ops. Numerical experiments which show the effectiveness of our algorithm are presented.

Keywords: Band Toeplitz systems, Newton’s iteration, regularization, displacement rank, image processing.

1. INTRODUCTION

Solving linear systems with noise is commonly encountered in image restoration problems. Typically, an image is represented by an m × n array X = (x_{i,j}) where the entry x_{i,j} provides the value of the pixel with coordinate (i, j) of the image. Using a vector notation, it is more convenient to represent X as a vector x = vec(X) obtained by arranging columnwise the entries of X. In this way, the typical problem of image restoration is to solve the linear system

Ax = b − η  (1)

where x is the original image, A is the blur operator so that Ax is the blurred image, and η is the additive noise of which we have a partial information, like some statistical properties or its norm. The vector b = vec(B) represents the blurred image B with the additive noise so that, solving the linear system

Ay = b  (2)

provides the vector y = x + A^{−1}η which may differ very much from x due to the noise and to the very large ill conditioning of A, that is typical for blur operators.

From the space invariance of the point spread function (PSF) that defines the matrix A and from its local action (a single source of light is blurred into a small spot wherever the source is located), the matrix A inherits the 2-level Toeplitz structure and the 2-level band structure. More precisely, if M = (m_{i,j})_{i,j=−k,k} denotes the mask representing the PSF, i.e., such that a unit pixel of coordinates (r, s) is spread onto the pixels m_{r+i,s+j}, i, j = −k, k, then the matrix A is a block Toeplitz matrix with a (2k +1)-diagonal block banded structure where the general block A_{i} on the i-th block diagonal, for i = −k, …, k, is a (2k + 1)-diagonal banded Toeplitz matrix with entry m_{i,j} on the j-th diagonal for j = −k, …, k.

In certain cases the matrix A can be written as T_{1} ⊗ T_{2}, where ⊗ denotes the Kronecker product. We refer to this situation as to the case of separable PSF. An important instance of separable PSF occurs if the entries of the mask are defined by means of a Gaussian function, say, m_{i,j} = β exp(−α(i^{2} + j^{2})) for given α, β > 0. In fact, T_{1} = T_{2} = (t_{j−i}), t_{h} = \sqrt{β} exp(−αh^{2}), h = −k, …, k, t_{h} = 0 elsewhere.

bini@dm.unipi.it, farusi@iei.pi.cnr.it, fiorent@dm.unipi.it, meini@dm.unipi.it

A fundamental problem in this context is to design algorithms for the solution of (2) providing a good approximation of the true image x, solution of (1), at a low computational cost. While the latter feature relies on exploiting the specific structures of the matrix A, the former is usually based on suitable regularization techniques (like Tikhonoff regularization), or on regularizing features of the conjugate gradient iterations (see Ref. 1).

Here, we introduce a new method for the regularized solution of the system (2) based on Newton’s iteration for matrix inversion and on the concept of approximate displacement rank of Ref. 2,3.

The new algorithm is particularly effective in the case of a separable PSF but can be used as a tool for speeding up the convergence, still keeping the regularizing properties, in the general case of nonseparable PSFs. The regularizing features are provided by the use of Newton’s iteration that is a well known quadratically convergent method for the approximation of the inverse matrix. In fact, the convergence of Newton’s iteration is much faster in the subspaces which do not contain the undesired noise components and is slow along the remaining components. This feature can be fine tuned and even amplified with the choice of a suitable initial approximation. We denote this choice as preconditioning the Newton iteration since its effect is somehow similar to the one of the preconditioned conjugate gradient.

The regularization property of Newton’s iteration seems to have been neglected in the current and past literature. Probably the reason is that a single iteration step is so expensive that the method is actually unpractical. In fact for an n × n matrix the cost of one step of Newton’s iteration amounts to O(n^{3}) arithmetic operations (ops). However, using the tools introduced in Ref. 2,3 and in particular the concept of approximate displacement rank, we are able to dramatically reduce the cost of a single step applied to a Toeplitz matrix to O(n log n) ops. This makes Newton’s iteration an effective tool for the design of regularizing algorithms.

Roughly speaking, our technique provides a displacement representation (see Ref. 4) of “regularized” inverses T_{1}^{+} and T_{2}^{+} of T_{1} and T_{2}, respectively. In this way, if we denote X and B the images corresponding to x and b (i.e., such that vec(X) = x and vec(B) = b, respectively), then a regularized approximation X̂ of X, solution of (1), is obtained as X̂ = T_{2}^{+}B(T_{1}^{+})^{T}. Moreover, the displacement structures of T_{1}^{+} and T_{2}^{+} allow one to compute the products T_{2}^{+}B(T_{1}^{+})^{T} with only O(n^{2} log n) ops, where for simplicity here and hereafter we assume m = n. The computation of T_{1}^{+} and T_{2}^{+} is performed in just O(n log n) ops.

For nonseparable problems, we use the same technique inside an iterative procedure for the successive approximation of the solution, based on the approach of Ref. 5,6 which relies on the singular values decomposition (SVD) of the mask M. More precisely, from the SVD M = ∑_{j=1}^{2k+1} σ_{j}u_{j}v_{j}^{T}, where σ_{1} ≥ … ≥ σ_{2k+1} ≥ 0 are the singular values and u_{j} = (u_{j,h})_{h=−k,k}, v_{j} = (v_{j,h})_{h=−k,k} are the left and right singular vectors, respectively, we obtain the splitting A = T_{1}^{(1)} ⊗T_{2}^{(1)} +∑_{j=2}^{2k+1} T_{1}^{(j)} ⊗T_{2}^{(j)}, where T_{1}^{(j)} = \sqrt{σ_{j}}U_{j}, T_{2}^{(j)} = \sqrt{σ_{j}}V_{j}, U_{j} and V_{j} are the (2k+1)-diagonal Toeplitz matrices defined by the vectors u_{j} and v_{j}, respectively, i.e., the entries on the i-th diagonal of U_{j} and V_{j} are u_{j,i} and v_{j,i}, respectively, for i = −k, …, k. For the PSFs encountered in practical problems the most part of information is contained in σ_{1}u_{1}v_{1}^{T} (see Ref. 5,6) whereas the singular values σ_{i}, for i > 1, are relatively negligible. This means that T_{1}^{(1)} ⊗ T_{2}^{(1)} is already a good approximation of A and can be used as a preconditioner in PCG iterations or in Jacobi-like iterations. Concerning the latter, the solution of the system is approximated with the iteration scheme

x_{i+1} = x_{i} + (T_{1}^{(1)}^{+} ⊗ T_{2}^{(1)}^{+}) (b − Ax_{i}) , i ≥ 0.  (3)

The paper is organized as follows. In section 2 we introduce the tools needed for the design of our algorithm including Newton’s iterations, regularization properties expressed in terms of filters, displacement rank and ϵ-displacement rank, orthogonal displacement representations. In section 3 we describe our modified version of Newton’s iteration and apply it to the separable and to the nonseparable case. Finally in section 4 we describe and comment our numerical experiments.

2. THE TOOLS

2.1. Newton’s Iteration

Newton’s iteration is a valid tool for solving problems whose solution can be formally viewed as the zero of a known function. This method has been analyzed in sections 1.6, 2.11, 3.5 of Ref. 7 and, concerning Toeplitz matrices has been investigated in Ref. 8,2,3.

Let A be a nonsingular matrix and consider the sequence of matrices generated, from an assigned X_{0}, in the following way

X_{i+1} = 2X_{i} − X_{i}AX_{i}.  (4)

This sequence 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.

Denoting R_{i} = I − X_{i}A the residual error at the i-th step, from (4) we obtain that R_{i+1} = R_{i}^{2}, that is, R_{i} = R_{0}^{2^{i}}. The above relation provides information on the convergence of the general iteration (4).

If A ∈ R^{n×n} is symmetric and positive definite with eigenvalues 0 < λ_{n} ≤ λ_{n−1} ≤ … ≤ λ_{1} = ρ(A), then the choice X_{0} = θI, with θ ≤ 1/ρ(A), say, θ = 1/||A|| where ||·|| is any operator norm, provides a sequence X_{i} converging to A^{−1}. In fact, the matrix R_{0} = I − θA is symmetric, nonnegative definite with eigenvalues 0 ≤ 1 − θλ_{1} ≤ … ≤ 1 − θλ_{n} < 1, therefore ρ(R_{0}) = 1 − θλ_{n} ≥ 1 − 1/cond_{2}(A), cond_{2}(A) = ||A||_{2}||A^{−1}||_{2} and ρ(R_{i}) = (1 − θλ_{n})^{2^{i}}. In this way, since (1 − 1/x)^{x} < 1/e < 1/2, after i = ⌈log_{2}(1/(θλ_{n}))⌉ iterations, we have ρ(R_{i}) ≤ 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 (4) for

ν = ⌈log_{2}(1/(θλ_{n}))⌉ + ⌈log_{2} log ϵ^{−1}⌉ ≥ ⌈log_{2} cond_{2}(A)⌉ + ⌈log_{2} log ϵ^{−1}⌉

steps.

The first additive term in the above expression, that in Ref. 2 is called overhead cost or overhead number of iterations, provides the number of steps needed in order to arrive at a pure double exponential convergence of order 2, where for pure double exponential convergence of order 2 we mean that the error tends to zero as a^{2^{i}} for a ≤ 1/2. The second term provides the number of steps of pure double exponential convergence that is needed to reach the residual error bound ϵ. A better choice of θ that still guarantees convergence is θ = 2/||A||.

For a generic matrix A, the starting approximation X_{0} = θA^{T} with θ ≥ ||A^{T} A||, guarantees the convergence. In fact, R_{0} = I − θA^{T} A is such that ρ(R_{0}) < 1. Moreover, if A = UΣV^{T} = ∑_{j=1}^{n} σ_{j}u_{j}v_{j}^{T} is the SVD of A, where Σ = diag(σ_{1}, …, σ_{n}), σ_{1} ≥ … ≥ σ_{n} ≥ 0, we have ρ(R_{0}) = 1 − θσ_{n}^{2} ≥ 1 − 1/cond_{2}(A)^{2}. Thus, the overhead cost becomes 2⌈log_{2} cond_{2}(A)⌉. That is, with respect to the positive definite case with the choice X_{0} = θI, the overhead cost is doubled. Also for a general matrix A we may choose θ = 2/||A^{T} A|| and convergence is guaranteed.

2.2. Regularizing Properties of Newton’s Iteration

Let Ax = b − η be the system where the noise η has been removed from the known term vector b, and let Ay = b the system with the noise. Indeed, y = x + A^{−1}η. From the SVD of A, A = UΣV^{T}, we easily deduce that

y = x + ∑_{j=1}^{n} \frac{u_{j}^{T} η}{σ_{j}} v_{j},  (5)

that is, the components of the error η along the vectors v_{i} are strongly amplified whenever σ_{j} ≈ 0. This amplification of the noise spoils the solution in its components along the correspondent vectors v_{j}.

Consider the approximation y_{i} = X_{i}b, where X_{i} is the i-th approximation to A^{−1} provided by Newton’s iteration with X_{0} = θA^{T}. Then we have

y_{i} = ∑_{j=1}^{n} ϕ_{j}^{(i)} \frac{u_{j}^{T} b}{σ_{j}}v_{j}  (6)

where the “filters” ϕ_{j}^{(i)} are such that

ϕ_{j}^{(i)} = 1 − (1 − θσ_{j}^{2})^{2^{i}}.  (7)

Moreover,

y_{i} = x − R_{i}x + (I − R_{i})A^{−1}η = x − ϵ_{i} + δ_{i},  (8)

where ϵ_{i} is the approximation error of Newton’s iteration and δ_{i} is the contribution of the noise η.

We easily find that

ϵ_{i} = R_{i}x = ∑_{j=1}^{n}(1 − θσ_{j}^{2})^{2^{i}}(v_{j}^{T} x)v_{j}, δ_{i} = ∑_{j=1}^{n} ϕ_{j}^{(i)} \frac{u_{j}^{T} η}{σ_{j}} v_{j}.  (9)

Figure 1. Functions ϕ^{(i)}(t) = 1 − (1 − t^{2})^{2^{i}} for i = 1, …, 5 over [0, 1]

In particular we point out that, for i → +∞, the convergence to zero of the approximation error ϵ_{i} is faster for the components along v_{j} such that σ_{j} ≈ 0 than for the ones corresponding to larger σ_{j}’s. In other words, the smaller is σ_{j}, the larger is the overhead cost along the j-th component, thus the pure double exponential convergence starts earlier for the larger σ_{j}’s.

Concerning the noise, we point out that the amplification along the component v_{j} given by the factor \frac{1}{σ_{j}} in (5),

is damped in (9) by the “filter” ϕ_{j}^{(i)} that is naturally introduced by Newton’s iteration. In order to better appreciate the role of this filter, let us observe that for the Euclidean norm || · ||_{2} it holds ||A^{T} A||_{2} = σ_{1}^{2}, thus for θ = \frac{1}{||A^{T} A||_{2}}

we have θσ_{j}^{2} = (\frac{σ_{j}}{σ_{1}} )^{2} < 1, moreover, if θ = \frac{2}{||A^{T} A||_{2}} then θσ_{j}^{2} = 2 (\frac{σ_{j}}{σ_{1}} )^{2}. Therefore we may plot the continuous

counterpart ϕ^{(i)}(t) = 1 − (1 − t^{2})^{2^{i}} of ϕ_{j}^{(i)} obtained by replacing \sqrt{θ}σ_{j} with the continuous variable t, where t ranges in [0, 1], or in [0, \sqrt{2}] if θ = 2/σ_{1}^{2}. From figure 1 we may see that the components of the noise corresponding to the values of t close to zero (i.e., the most dangerous ones) are strongly filtered since ϕ^{(i)}(t) is close to zero, whereas the remaining ones are almost left unchanged.

We may estimate the threshold value under which the filter acts. We find that for t < t̂(i) = \sqrt{1 − \frac{1}{2^{2^{i−1}}}} the filter takes values less than 1/2. We may choose the number of steps of Newton’s iteration according to this threshold. Another possible criterion for determining the number of steps is based on the following observation. From (9) it follows that the noise is amplified by the factor ψ^{(i)}(t) = ϕ^{(i)}(t)/t, moreover, if i < log_{2}(1/t̃) then we have ψ^{(i)}(t) < 1 for 0 ≤ t ≤ t̃. That is, after ν steps of Newton’s iteration, the noise is not amplified along the components for which σ_{i} < 1/2^{ν}. Once we have chosen the cut-off value t̃, based on suitable information about the noise, we may determine the number of steps of Newton’s iteration which guarantee that the noise is not amplified under the cut-off level.

It is interesting to observe that the expressions (6), (7), (8) and (9) coincide with the analogous of the Landweber iteration (see Ref. 1, page 137) with the difference that the exponent in the filters ϕ_{j}^{(i)} in the Landweber iteration is just i, whereas in the Newton iteration is 2^{i}. This reflects the quadratic convergence of Newton’s iteration versus the linear convergence of Landweber’s iteration. Similar regularization properties are obtained with X_{0} = θI.

It is possible to improve the convergence of the Newton iteration still keeping the regularization properties, by choosing a different initial approximation X_{0}. Let us consider the case where A is positive definite and observe that if C is positive definite then C^{1/2} is well defined and the matrix CA is similar to the symmetric positive definite matrix C^{1/2}AC^{1/2}. Therefore, by choosing X_{0} = θC with θ ≤ 1/||CA||, the residual R_{0} = I − X_{0}A = I − θCA has eigenvalues 1 − θλ_{j} where 0 ≤ λ_{1} ≤ … ≤ λ_{n} are the eigenvalues of CA, and the same arguments used in the analysis of the case X_{0} = θI apply. In fact, it holds

R_{i} = ∑_{j=1}^{n}(1 − θλ_{j})^{2^{i}}u_{j}v_{j}^{T} ,  (10)

where u_{j}, v_{j} are the right and left eigenvectors of CA, respectively.

Similarly, for a general matrix A we may choose X_{0} = θA^{T} C^{T} C, for any matrix C, θ ≤ 1/||CA||^{2}. In this case we have R_{0} = I − θ(CA)^{T} (CA) and ρ(R_{0}) < 1. Alternative choices for X_{0} are X_{0} = θCC^{T} A^{T} or X_{0} = \frac{θ}{2}(CC^{T} A^{T} + A^{T} C^{T} C). A similar analysis can be carried out and the overhead cost now becomes 2 log cond_{2}(CA). Thus the convergence can be improved if cond_{2}(CA) << cond_{2}(A). The analogy with the preconditioned conjugate gradient technique where the rate of (linear) convergence depends on cond_{2}(CA) (for the preconditioned version) and on cond_{2}(A) (for the unpreconditioned version), leads us to define the matrix C as the preconditioner and to refer to Newton’s iteration with X_{0} depending on C as to preconditioned Newton’s iteration.

The wide literature on the design of preconditioners (see, for instance Ref. 9) provides us many techniques for the construction of a suitable C. However, the choice of C is rather delicate since it may strongly accelerate the convergence but, at the same time, it may destroy the regularization features of Newton’s iteration. In particular, this holds for the circulant preconditioner of R. Chan^{9} unless its inversion is replaced by some sort of generalized inversion.

For a general n × n matrix A, Newton’s iteration has the cost of O(n^{3}) ops per step. This is the main drawback that makes this iteration unpracticable for the effective inversion of a matrix. However, for special classes of matrices, say, for Toeplitz matrices, the cost can be reduced to O(n log n) ops per step. This can be accomplished by using the concept of approximate displacement rank.

2.3. Approximate Displacement Rank

The concept of displacement operators and displacement rank, introduced in Ref. 4,10, and successively elaborated by other authors (see Ref. 11 and section 2.11 of Ref. 7), is a powerful tool for dealing with Toeplitz matrices. In certain situations this tool seems to be inadequate to cope with the specific computational issues encountered in practice and a modified version of it with better “numerical” features seems best suited for several applications. Here we recall the main results concerning displacement rank, orthogonal displacement representation and approximate displacement rank and prepare the tools for the effective implementation of Newton’s iteration. Throughout this section we will follow Ref. 2,3,12.

Let us recall the operator

Δ(A) = AZ − ZA, Z = (z_{i,j}), z_{i+1,i} = 1, z_{i,j} = 0, otherwise,  (11)

applied to n × n matrices A, and define the displacement rank of A as drk(A) = rank(Δ(A)). A Toeplitz matrix has displacement rank at most 2, moreover, drk(A) is a sort of measure of the distance of A from the algebra of lower triangular Toeplitz matrices, for which the displacement rank is zero.

Define L(a) = ∑_{i=0}^{n−1} a_{i+1}Z^{i} the n × n lower triangular Toeplitz matrix having a = (a_{i})_{i=1,n} as first column. It holds (see section 2.11 of Ref. 7)

Δ(A) = ∑_{i=1}^{r} u_{i}v_{i}^{T} ⇒ A = L(Ae_{1}) + ∑_{i=1}^{r} L(u_{i})L^{T} (Zv_{i}),  (12)

where u_{i}, v_{i}, i = 1, …, r, are n-dimensional vectors, e_{1} denotes the first column of the identity matrix. The expression in the right-hand side of (12) is called the displacement representation of A and the vectors a = Ae_{1}, u_{i}, v_{i}, i = 1, …, r are called a set of generators for the displacement representation of A. A matrix A with displacement rank r is called Toeplitz-like if r << n.

A simple consequence of the relation

Δ(A^{−1}) = −A^{−1}Δ(A)A^{−1},  (13)

that can be directly derived from (11), is that drk(A^{−1}) = drk(A) and the following displacement representation of A^{−1} holds

A^{−1} = L(A^{−1}e_{1}) − ∑_{i=1}^{r} L(A^{−1}u_{i})L^{T} (ZA^{−T} v_{i}).  (14)

Since the product of an n × n Toeplitz matrix and a vector can be performed in O(n log n) ops by means of FFT, the representations (12), (14) allow us to make elementary operations among Toeplitz-like matrices at a low cost. In particular, computing the product of a matrix with displacement rank r and a vector costs just O(rn log n) ops.

The set of generators for the displacement representation of A is not unique. Indeed, observe that, if Δ(A) = UV^{T}, for U, V ∈ R^{n×r} then Δ(A) = ÛV̂^{T} for Û = US, V̂ = V S^{−T}, where S ∈ R^{r×r} is any nonsingular matrix. Therefore another displacement representation can be derived by the latter decomposition. If the product of the norms of the columns of U and the corresponding columns of V is much larger than the norm of A, then computing the matrix vector product Ax by means of the displacement representation formula (12) may lead to large numerical cancellation errors. This drawback can be overcome by using the concept of orthogonal displacement representation.^{2,3}

Definition 2.1. If Δ(A) = UΣV^{T}, Σ = diag(σ), σ = (σ_{1}, …, σ_{r}), U, V ∈ R^{n×r}, is the singular values decomposition of Δ(A), then the quadruple (a, U, σ, V) defines a set of generators, such that

A = L(a) + ∑_{i=1}^{r} σ_{i}L(u_{i})L^{T} (Zv_{i}),  (15)

where u_{i} = Ue_{i}, v_{i} = V e_{i}, i = 1, …, r, and e_{i} denotes the i-th column of the n × n identity matrix. We refer to this particular set of generators as to an orthogonal displacement generator for the displacement representation of A, and we set ODG(A) = (a, U, σ, V). The matrix representation (15) is called orthogonal displacement representation of A and we write A = ODR(a, U, σ, V).

Observe that for the Euclidean norm || · || we have ||Δ(A)|| = ||AZ − ZA|| ≤ 2||A||, whence for an orthogonal set of generators we have ||U|| = ||V || = 1, σ_{1} = ||Δ(A)|| ≤ 2||A||. This fact overcomes the possibility of numerical cancellation.

A matrix A which is symmetric with respect to the secondary diagonal, i.e., a_{i,j} = a_{n−j+1,n−i+1}, is called persymmetric. A persymmetric matrix A is such that A^{T} = JAJ, where J is the permutation matrix having 1 in the anti-diagonal. An interesting example of persymmetric matrices is the class of Toeplitz matrices. It can be easily proved that Δ(A)J is skew-symmetric if A is persymmetric. The following result holds:

Theorem 2.2. [Displacement representation of a persymmetric matrix]. If the persymmetric matrix A has displacement rank r, then r = 2h and there exist matrices W ∈ R^{n×r}, D = diag(d_{1}, …, d_{h}) ∈ R^{h×h} such that W^{T} W = I_{r}, |d_{1}| ≥ |d_{2}| ≥ … ≥ |d_{h}|, |d_{j}| = σ_{2j−1} = σ_{2j}, j = 1, …, h, where σ_{1}, …, σ_{r} are the singular values of Δ(A), 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}))  (16)

where a = Ae_{1}, w_{i} = W e_{i}, i = 1, …, r. An orthogonal generator for A is completely defined by the triple (a, W, d) where d = (d_{i})_{i,1,...,h}.

This fact can be exploited for reducing the computational cost of multiplying the matrix A and a vector. Indeed, it is easy to prove that, given the FFTs associated with x and y, only the FFTs associated with w, x and y are sufficient for the computation of the products L(w)x and L(ZJw)y.

Definition 2.3. If A is persymmetric and Δ(A) = W (D ⊗ S)W^{T} J, we define the triple (a, W, d), where a = Ae_{1}, d = (d_{i})_{i=1,h} a persymmetric orthogonal displacement generator of A and we use the same notation ODG(A) = (a, W, d). Similarly, the representation

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

is called persymmetric orthogonal displacement representation of A and we use the same notation A = ODR(a, W, d).

The following algorithms provide a detailed description of the computation of the ODG of a matrix.

Algorithm 2.4. [Computing an ODG of a persymmetric matrix]. Input: Matrices U ∈ R^{n×k}, M ∈ R^{k×k} such that Δ(A) = F MF^{T} J, M = −M^{T}, where A is persymmetric; Output: ODG(A) = (a, W, d). Computation:

– Compute the QR factorization F = QR.

– Compute H = RMR^{T}.

– Compute the SVD H = UΣV.

– Compute the sign γ_{j} = sign(u_{i,2j−1}, v_{i,2j}), j = 1, …, h, where i is such that u_{i,2j−1}v_{i,2j} ≠ 0.

– Set W = QU, d_{i} = γ_{i}σ_{2i−1}, i = 1, …, h, d = (d_{i}), a = Ae_{1}.

Now we recall the concept of approximate displacement rank, or ϵ-displacement rank, introduced and analyzed in Ref. 2,3.

Definition 2.5. For a given ϵ > 0 and for the Euclidean norm || · ||, define the ϵ-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). If r = drk_{ϵ}(A) (r = rdrk_{ϵ}(A)), and if Δ(A) = ∑_{i=1}^{r} u_{i}v_{i}^{T} + E, with ||E|| ≤ ϵ (||E|| ≤ ϵ||A||), we say that the matrix

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

is an ϵ-displacement (relative ϵ-displacement) representation of A.

The following properties hold (see Ref. 3,12).

Theorem 2.6. Let Δ(A) = UΣV^{H} be the singular value decomposition of Δ(A), where σ_{1} ≥ … ≥ σ_{k} > 0 are the singular values of Δ(A). Let σ_{1} > ϵ ≥ σ_{k}; then drk_{ϵ}(A) = r if and only if σ_{r} > ϵ ≥ σ_{r+1}. Similarly, let σ_{1} > ϵ||A|| ≥ σ_{k}; then rdrk_{ϵ}(A) = r if and only if σ_{r} > ϵ||A|| ≥ σ_{r+1}.

In the hypotheses of Theorem 2.6 we call the quadruple (a, U_{ϵ}, σ_{ϵ}, V_{ϵ}) orthogonal (relative) ϵ-displacement generator of A if σ_{ϵ} = (σ_{1}, …, σ_{r}), U_{ϵ} and V_{ϵ} are made up by the first r columns of U and V, respectively, and a = Ae_{1}. Moreover, the representation

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

is called orthogonal (relative) ϵ-displacement representation of A.

We refer the reader to Ref. 2,3,12, for properties and results concerning the approximate displacement rank.

In order to construct orthogonal ϵ-displacement generators for a matrix 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̂), Û = [Ue_{1}, …, Uek̂], V̂ = [V e_{1}, …, V ek̂] ∈ R^{n×k̂} \\ σk̂ > σ_{1}ϵ ≥ σ_{k̂+1}, \end{matrix}  (17)

for a general A, and in the following way for a persymmetric matrix A

\begin{matrix}Trunc_{ϵ}(a, W, d) = (a, Ŵ, d̂), \\ d̂ = (d_{1}, …, d_{k̂/2}), Ŵ = [W e_{1}, …, W ek̂] ∈ R^{n×k̂} \\ |dk̂| > |d_{1}|ϵ ≥ |d_{k̂+1}|. \end{matrix}  (18)

Equations (17), (18) and Algorithm 2.4 provide a means for computing an ϵ-ODG of A whenever any orthogonal displacement representation is available.

3. NEWTON’S ITERATION FOR TOEPLITZ MATRICES

With the tools introduced in the previous section we are able to define our regularized version of Newton’s iteration adjusted to cope with the case of Toeplitz matrices, based on Ref. 2,3. Local convergence analysis of a similar modification of Newton’s iteration has been performed in Ref. 8.

If X_{0} is such that X_{0} = JX_{0}^{T} J, then all the matrices X_{i} satisfy the equation X_{i} = JX_{i}^{T} J and we may represent each X_{i} by means of its persymmetric orthogonal displacement generator of Definition 2.3. Observe that the displacement rank of X_{i} may grow as i tends to infinity, even though the limit A^{−1} of X_{i} has displacement rank 2. However we may overcome this drawback by means of the ϵ-displacement rank. In fact, for any ϵ > 0, the value drk_{ϵ}(X_{i}) is asymptotically bounded by drk(A^{−1}).

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}  (19)

where Trunc_{ϵ}(·) is defined in (18) 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 ϵ_{i} that guarantees the convergence of the modified Newton method can be easily designed. We recall the following result of Ref. 2,3.

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

For θ = 1/cond_{2}(A), the bound on ϵ turns into ϵ ≤ ||R_{i}||^{2}/(4n(2h_{i}+k)σ_{1}cond_{2}(A)). Indeed, this bound guarantees the convergence whatever the matrix A is. However, due to the regularization properties, we do not need to perform many steps of Newton’s iteration, rather, the value of log_{2} θ^{−1} is customarily close to 1 so that our regularized version of Newton’s iteration stops at the very early stages.

In order to analyze the iteration (19) we will relate ODG(X) = (x, W_{X}, d_{X}) with ODG(Y) = (y, W_{Y} , d_{Y}), where Y = 2X − XAX.

It is immediate to verify that:

\begin{matrix}y = 2x − XAx \\ Δ(Y) = F MF^{T} J \\ M = [ \begin{matrix}2D_{X} & O & −D_{X} \\ 0 & −D_{A} & O \\ −D_{X} & O & O \end{matrix} ] ⊗ [ \begin{matrix}0 & 1 \\ −1 & 0 \end{matrix} ] \\ F = [ W_{X} XW_{A} XAW_{X} ] \end{matrix}  (20)

Recovering ODG(Y) from (20) is almost inexpensive if performed by means of Algorithm 2.4. Now we are ready to describe our implementation of Newton’s iteration based on (19), where we use the infinite norm of a matrix A defined by ||A||_{∞} = max_{i} ∑_{j} |a_{i,j}|.

Algorithm 3.2. [Regularized Newton’s iteration with ϵ-displacement rank].

Input: Integer n, the n × n Toeplitz matrix assigned in terms of its first row and its first column; a strategy for the choice of the truncation value ϵ_{i}; the number k of steps.

Output: A regularized approximation X of A^{−1} given in terms of its ODG (x, W_{X}, d_{X}).

Computation:

– Compute θ = 1/||A||_{∞}||A^{T} ||_{∞}, set X_{0} = θA^{T}, i = 0 and compute ODG(A) and ODG(A^{T}).

– For i = 0, 1, …, k − 1 compute (y, W_{Y} , d_{Y}) = ODG(Y), where Y = 2X_{i} − X_{i}AX_{i}, by means of (20) and Algorithm 2.4.

– Compute ODG(X_{i+1}) = Trunc_{ϵ_{i}}(ODG(Y)), by means of (18).

– output ODG(X_{k}).

The algorithm can be easily modified by introducing the preconditioning stage where X_{0} = θCC^{T} A^{T} for a suitable matrix C having low displacement rank, or X_{0} = \frac{θ}{2}(CC^{T} A^{T} + A^{T} C^{T} C) if we want to preserve the persymmetry of the iterates.

3.1. The Separable Case

Consider the system (2) where A = T_{1} ⊗ T_{2} and T_{1}, T_{2} are (banded) Toeplitz matrices. This case is encountered in practice, for instance, for the PSF given by the mask M = (m_{i,j})_{i,j=−k,k}, where m_{i,j} = 1/α, α = (2k + 1)^{2}; in this case T_{1} = T_{2} is a 2k + 1 diagonal Toeplitz matrix whose nonzero entries are equal to 1/(2k + 1). Another interesting case is given by the mask M = (m_{i,j}) defined by a Gaussian function, i.e., m_{i,j} = β exp(−αi^{2} − αj^{2}). In fact, in this case T_{1} = T_{2} is the (2k + 1)-diagonal Toeplitz matrix whose nonzero entries are \sqrt{β} exp(−αi^{2}).

By representing the unknown vector y by means of a matrix Y and b by B where y = vec(Y), b = vec(B), the system Ay = b turns into T_{2}Y T_{1}^{T} = B, from which we have Y = T_{2}^{−1}BT_{1}^{−T}. A regularized approximation of Y is given by T_{2}^{+}BT_{1}^{+}, where T_{2}^{+} and T_{1}^{+} are regularized approximations to T_{1}^{−1} and T_{2}^{−1}, respectively, obtained by applying a suitable (small) number of steps of Newton’s iteration to T_{1} and T_{2}. The computation of T_{2}^{+}BT_{1}^{+} consists in performing 2 products between a Toeplitz-like matrix and a general n × n matrix (if the image has size n × n); this computation requires O(n^{2} log n) arithmetic operations.

3.2. The General Case

Let M = ∑_{j=1}^{2k+1} σ_{j}u_{j}v_{j}^{T} be the SVD of the mask M associated with the PSF matrix A. Let T_{1}^{(j)} = \sqrt{σ_{j}}U_{j},

T_{2}^{(j)} = \sqrt{σ_{j}}V_{j}, where U_{j} and V_{j} are the (2k + 1)-diagonal Toeplitz matrices defined by the vectors u_{j} and v_{j}, respectively. Let (T_{1}^{(1)})^{+} and (T_{2}^{(1)})^{+} be the regularized approximations to T_{1}^{(1)} and T_{2}^{(1)} provided by a suitable number of steps of Newton’s iteration. Consider the splitting

A = ((T_{1}^{(1)})^{+} ⊗ (T_{2}^{(1)})^{+})^{−1} − (((T_{1}^{(1)})^{+} ⊗ (T_{2}^{(1)})^{+})^{−1} − A)

and observe that the iterative scheme associated with this partitioning is (3); moreover, the iteration matrix is P = I − ((T_{1}^{(1)})^{+} ⊗ (T_{2}^{(1)})^{+}) A. If ρ (((T_{1}^{(1)})^{+} ⊗ (T_{2}^{(1)})^{+})A) < 1 the iterative scheme (3) converges. If the singular values and the singular vectors of A do not differ much from the corresponding singular values and the singular vectors of T_{1}^{(1)} ⊗ T_{2}^{(1)} then the iterative scheme has regularizing properties.

4. NUMERICAL EXPERIMENTS

We have applied our technique to the case where A is the two-level banded Toeplitz matrix defined by the 31 × 31 mask m_{i,j} = βe^{−α(i^{2}+j^{2})}, i, j = −15, 15, with α = 0.1, and β so that ∑_{i,j} m_{i,j} = 1. In this way it holds A = T ⊗ T, T = (t_{i,j}), t_{i,j} = \sqrt{β}e^{−α(i−j)^{2}} if |i − j| ≤ 15, t_{i,j} = 0 elsewhere. We set n = m = 512, i.e., our image is formed by 512^{2} pixels, where each pixel is represented by an integer between 0 and 255. We have considered the case where the noise η is originated by the numerical rounding to the integer part of the blurred image (numerical noise) which has real entries.

We have applied Newton’s iteration with the parameter θ = 1, in fact it holds σ_{1} ≤ 1, and with the initial approximation X_{0} = θI and X_{0} = T^{T}. The results of the numerical experiments are reported in figure 2.

The regularized features of Newton’s iteration are evident from figure 2: in fact, three iterations with X_{0} = I are enough to provide a good approximation of the image still keeping the noise under a suitable level, while a larger number of iterations amplifies the error. With the choice X_{0} = T^{T} the number of iterations sufficient to recover the image is larger (according to our theoretical estimates): 8 iterations are enough to recover the image without much amplifying the noise. With 10 iterations the noise slightly appears in the image. It is interesting that the

Original image  Blurred image

3 Newton’s iterations with X_{0} = I  4 Newton’s iterations, X_{0} = I

8 Newton’s iterations, X_{0} = T^{T}  10 Newton’s iterations, X_{0} = T^{T}

Figure 2. Separable PSF

Blurred image  it = 1, k = 3, θ = 0.1

it = 2, k = 3, θ = 0.1  it = 4, k = 3, θ = 0.1

it = 1, k = 6, θ = 0.1  it = 2, k = 6, θ = 0.1

Figure 3. Nonseparable PSF amplification of the noise in the cases X_{0} = I and X_{0} = T^{T} occurs with different features, in particular the noise is coarse-grained if X_{0} = T^{T}.

We have applied our technique to the reconstruction of the image of Saturn blurred with a nonseparable PSF, that can be obtained by ftp from ftp.stsci.edu (see also Ref. 13,1).

We report the images obtained by applying it iteration of our scheme (3) where the regularized inversion of the matrices T_{1}^{(1)} and T_{2}^{(1)} is obtained by means of k Newton’s iterations with the initial point X_{0} = θI. The results are reported in figure 3. We have chosen θ = 0.1 since a larger value, say θ = 1 still guaranteeing convergence, would have caused a too fast convergence of Newton’s iteration with a consequent amplification of the noise even with a single step. As we can see from figure 3, the image obtained after 4 iterations of (3), where the regularized inverses of T_{1}^{(1)} and T_{2}^{(1)} are obtained with 3 Newton’s steps and with X_{0} = θI, θ = 0.1, is the one which provides a good reconstruction with a very small noise. If we use more accurate approximations of the inverses of T_{1}^{(1)} and T_{2}^{(1)}, say with k = 6 Newton’s steps, the amplification of the noise is very strong.

REFERENCES

1. M. Bertero and P. Boccacci, Introduction to inverse problems in imaging, Institute of Physics Publishing, Bristol, 1998.

2. D. A. Bini and B. Meini, “Solving block banded block Toeplitz systems with structured blocks: algorithms and applications,” in Structured matrices: recent developments in theory and computation, D. A. Bini, E. Tyrtyshnikov, and P. Yalamov, eds., NOVA, 2000. In press.

3. D. A. Bini and B. Meini, “Approximate displacement rank and applications,” in Proceedings AMS Conference ”Structured Matrices in Operator Theory, Numerical Analysis, Control, Signal and Image Processing”, Boulder, V. Olshevsky, ed., AMS, 1999. to appear.

4. T. Kailath, S. Kung, and M. Morf, “Displacement ranks of matrices and linear equations,” J. Math. Anal. Appl. 68, pp. 395–407, 1979.

5. J. Kamm and J. G. Nagy, “Kronecker product and SVD approximations in image restoration,” Linear Algebra Appl. 284(1-3), pp. 177–192, 1998. ILAS Symposium on Fast Algorithms for Control, Signals and Image Processing (Winnipeg, MB, 1997).

6. J. Kamm and J. G. Nagy, “Optimal Kronecker product approximation of block Toeplitz matrices,” SIAM J. Matrix Anal. Appl. 22(1), pp. 155–172, 2000.

7. D. A. Bini and V. Pan, Matrix and Polynomial Computations, Vol. 1: Fundamental Algorithms, Birkhäuser, Boston, 1994.

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

9. R. H. Chan and M. K. Ng, “Conjugate gradient methods for Toeplitz systems,” SIAM Review 38, pp. 427–482, 1996.

10. T. Kailath, S. Kung, and M. Morf, “Displacement rank of a matrix,” Bulletin of the American Mathematical Society 1, pp. 769–773, 1979.

11. G. Heinig and K. Rost, Algebraic Methods for Toeplitz-like Matrices and Operators, Akademie-Verlag, Berlin, and Birkhäuser, Boston, 1984.

12. D. A. Bini and B. Meini, “Solving block banded block Toeplitz systems with banded Toeplitz blocks,” in Advanced Signal Processing Algorithms, Architectures, and Implementations IX, F. T. Luk, ed., vol. 3807, pp. 300–311, Proceedings of SPIE, 1999.

13. J. G. Nagy, “Decomposition of block Toeplitz matrices into a sum of Kronecker products with applications in image restoration,” Tech. Rep. 96-1, SMU Department of Mathematics, 1996.