Pivoted Cauchy-Like Preconditioners for Regularized Solution of Ill-Posed Problems.
Many ill-posed problems are solved using a discretization that results in a least squares problem or a linear system involving a Toeplitz matrix. The exact solution to such problems is often hopelessly contaminated by noise, since the discretized problem is quite ill conditioned, and noise components in the approximate null-space dominate the solution vector. Therefore we seek an approximate solution that does not have large components in these directions. We use a preconditioned conjugate gradient algorithm to compute such a regularized solution. A unitary change of coordinates transforms the Toeplitz matrix to a Cauchy-like matrix, and we choose our preconditioner to be a low rank Cauchy-like matrix determined in the course of Gu's fast modified complete pivoting algorithm. We show that  if the kernel of the ill-posed problem is smooth, then this preconditioner has desirable properties: the largest singular values of the preconditioned matrix are clustered around one, the smallest singular values, corresponding to the lower subspace, remain small, and the upper and lower spaces are relatively unmixed. The preconditioned algorithm costs only O(n lg n) operations per iteration for a problem with n variables.  The effectiveness of the preconditioner for filtering noise is demonstrated on three examples.
Introduction
. In fields such as seismography, tomography, and signal pro-
cessing, the process describing the acquisition of data can often be described by an
integral equation of the first kind
Z fi up
lo
where t denotes the kernel, -
f the unknown input function, and - g the output. When
it is appropriately discretized, the equation becomes a system of n linear equations of
the form
g:
In applications, properties of the kernel and the discretization process often cause T
to have a Toeplitz structure; that is,
constant along diagonals.
The discrete inverse problem is to recover -
f , given - g and T . However, the continuous
problem is generally ill-posed: i.e. small changes in - g cause arbitrarily large
changes in -
f . This is reflected in the discrete problem by ill-conditioning in the matrix
T . The recovery of -
f then becomes a delicate matter since the recorded data will
This work was supported by the National Science Foundation under Grant CCR 95-03126.
y Applied Mathematics Program, University of Maryland, College Park, MD 20742
z Department of Computer Science and Institute for Advanced Computer Studies, University of
Maryland, College Park, MD 20742 (oleary@cs.umd.edu).
likely have been contaminated by noise e. In this case, we have measured g rather
than - g, where
(1)
Due to the ill-conditioning of T and the presence of noise, exact solution of the linear
system will not lead to a reasonable approximation of -
f . Rather, regularization is
needed in order to compute an approximate solution f . Regularization can be thought
of as exchanging the original, ill-posed problem for a more well-posed problem whose
solution approximates the true solution. Many regularization methods, both direct
and iterative, have been discussed in the literature; see, for example, [12, 15, 9, 5]. In
this paper we will primarily be concerned with regularization via conjugate gradient
iterations [7, 22, 29], where the regularization parameter is the number of iterations.
Toeplitz matrices have several properties convenient for iterative methods like
conjugate gradients: multiplication of a Toeplitz matrix times a vector can be done in
O(n lg n) operations, and circulant preconditioners can be quite efficient [25, 3]. There
are some difficulties, though. The inverse of a Toeplitz matrix does not generally have
Toeplitz structure, and the fast factorization algorithms for Toeplitz matrices can
require as much as O(n 3 ) flops if pivoting is used to improve stability; see [27, 11, 4],
for example.
To overcome these difficulties, we make use of the fact that Toeplitz matrices are
related to Cauchy-like matrices by fast orthogonal transformations [17, 8, 10]. Cauchy-
like matrices, discussed in detail in x2, permit fast matrix-vector multiplication. But,
in contrast to Toeplitz matrices, the inverse of a Cauchy-like matrix is Cauchy-like,
and complete pivoting can be incorporated in its LDU factorization at a total cost of
The focus of this paper is the development of a Cauchy-like preconditioner that
can be used to accelerate convergence of the conjugate gradient iteration to a filtered
approximate solution of a problem involving a Toeplitz matrix. The regularizing
properties of conjugate gradients and our choice of preconditioner are discussed in x3.
Each iteration of our algorithm takes O(n lg n) operations, and computational issues
are discussed in x4. Section 5 contains numerical results and x6 presents conclusions
and future work.
2. Transformation from Toeplitz to Cauchy-like structure. A Cauchy-
like, or generalized Cauchy, matrix C has the form
' a T
1-i;j-n
(2)
It can also be defined as the unique solution of the displacement equation
where
a
The pair (A; B) is the generator of C with respect
to\Omega and \Theta, and ' - n is called
the displacement rank. For the matrices and displacement equations of interest here,
We exploit three important properties of Cauchy-like matrices.
Property 1. Row and column permutations of Cauchy-like matrices are Cauchy-
like, as are leading principal submatrices.
This property allows pivoting in fast algorithms for factoring Cauchy-like matrices
[17, 8].
Property 2. The inverse of a Cauchy-like matrix is Cauchy-like:
1-i;j-n
Heinig [17] gives an O(n lg 2 n) algorithm to compute X (with rows x T
(with rows w T
and\Omega\Gamma and explains how, using the FFT, a system
involving a Cauchy-like matrix can be solved in O(n lg 2 n). However, the algorithm
is very fragile. It can be unstable for large values of n and, even when used on a well
conditioned matrix, may require pivoting to maintain stability [18, 1]. Alternatively,
X and W can be determined from the relations
The third important property is that Toeplitz matrices also satisfy certain displacement
equations [21, 8] which allow them to be transformed via fast Fourier
transforms into Cauchy-like matrices [17, 8]:
Property 3. Every Toeplitz matrix T satisfies an equation of the form
R
The Toeplitz matrix T is orthogonally related to a Cauchy-like matrix
that satisfies the displacement equation
where
and F is the normalized inverse discrete Fourier transform matrix defined by
exp
(j
'-
1-j;k-n
Gohberg, Kailath, and Olshevsky [8] suggest a stable O('n 2 ) partial pivoting
algorithm to factor . Sweet and Brent [26] show, however, that element
growth in this algorithm depends not only on the magnitude of L and U , but on the
generator for the Cauchy-like matrix. For our test matrices, partial pivoting alone
did not provide the rank revealing information that we need.
Gu [10] presents an algorithm that can perform a fast O('n 2 ) variation of LU
decomposition with complete pivoting. Recall that in complete pivoting, at every
elimination step one chooses the largest element in the current submatrix as the
pivot in order to reduce element growth. Gu proposes instead that one find an entry
sufficiently large in magnitude by considering the largest 2-norm column of -
corresponding to the part that remains to be factored at each step. This algorithm
computes the factorization Alg. 2] using only the readily determined
generators (see x4), and Gu shows that it is efficient and numerically stable, provided
that element growth in the computed factorization is not large. For our purposes it
was convenient to set to obtain the equivalent
factorization
3. Regularization and preconditioning. If we wanted to solve the linear system
exactly, we would be finished: using the transformation to Cauchy-like
form and the fast factorization algorithms described above, computing this solution
would be an easy task. But the solution we seek is an approximate one, having noise
filtering properties, so we choose to use an iterative method called CGLS which, in
conjunction with an appropriate preconditioner, produces suitably filtered solutions.
Three assumptions will guide our analysis:
1. The matrix T has been normalized so that its largest singular value is of order
1.
2. The uncontaminated data vector -
satisfies the discrete Picard condition; i.e.,
the spectral coefficients of -
decay in absolute value like the singular values
[30, 14].
3. The additive noise is zero-mean white Gaussian. In this case, the components
of the error e are independent random variables normally distributed with
mean zero and variance ffl 2 .
We need to define the signal and noise subspaces. Using (1), let
V T be
the singular value decomposition of T , and expand the data and the noise in the basis
created by the columns of -
with
e. Under the white noise assumption, the coefficients j i
are roughly constant in size, while the discrete Picard condition tells us that the - fl i
go to zero at least as fast as the singular values oe i . Thus, components for which - fl i
is of the same order as j i are obscured by noise. Let m be such that j-fl i j AE jj i j for
we say that the
last
m columns of -
span the noise subspace, while the other columns span the
signal subspace. The basis for the signal subspace is further partitioned into the first
columns and the remaining -
which correspond to a transition subspace that
is generally difficult to resolve unless there is a gap in the singular value spectrum.
3.1. Regularization by preconditioned conjugate gradients. The standard
conjugate gradient (CG) method [19] is an iterative method for solving systems
of linear equations for which the matrix is symmetric positive definite. If the matrix
is not symmetric positive definite, one can use a variant of standard CG which
solves the normal equations in factored form. We refer to the resulting algorithm as
CGLS [19]. If the discrete Picard condition holds, then CGLS acts as an iterative
regularization method with the iteration index taking the role of the regularization
parameter [7, 13, 15]. Convergence is governed by the spread and clustering of the
singular values [28]. Therefore, preconditioning is often applied in an effort to cluster
the singular values, thus speeding convergence.
In the context of an ill-conditioned matrix T , we require a preconditioner for
CGLS which clusters the largest m singular values while leaving the small singular
values, and with them, the noise subspace, relatively unchanged. In this case, the first
few iterations of CGLS will quickly capture the solution lying within the subspace
spanned by the first m columns of V . A modest number of subsequent iterations will
provide improvement over the transition subspace, without significant contamination
from the noise subspace.
3.2. The preconditioner. Given the Toeplitz matrix T , let ~
its corresponding Cauchy-like matrix. Solving T is then equivalent to solving
~
Note that since F and S 0 are unitary matrices, then
~
that is,
C have the same singular values, and there is no mixing of
signal and noise subspaces.
A factorization of ~
C using a modified complete pivoting strategy may lead to an
interchange of rows (specified by a permutation matrix P ) and columns (specified by
a permutation matrix Q). Setting
the
problem we wish to solve is
We choose a preconditioner M for the left so that
and apply CGLS to the corresponding normal equations
Our choice of preconditioner M is derived from the leading m \Theta m submatrix of
Gu's modified complete pivoting LDU factorization of the matrix C as follows. Let
and write this equation in block form, where the upper left blocks are
m \Theta m:
Here are lower triangular, U 1 ; U 3 are upper triangular, and D 1 and D 2 are
diagonal. We choose as our preconditioner the matrix
I
I
I
3.3. Properties of the preconditioner. We begin with some theorems about
the clustering of the singular values of M \Gamma1 C. It is useful to decompose the matrix
into the matrix sum
using the block partitioning of the previous section.
be the sum of the absolute values of the entries in row i of C
be the largest of these quantities, and let - s be the largest such row sum for E 2 . The
case of interest to us is when these quantities are reasonably small.
We denote the k-th largest singular value of a matrix Z by oe k (Z), and the k-th
largest eigenvalue by - k (Z).
Theorem 3.1. If kC \Gamma1
largest singular values of M \Gamma1 C lie
in the interval
\Theta

Proof: The upper bound can be obtained by applying Gershgorin's theorem
[24][IV.2.1] to bound the eigenvalues of the matrix (M \Gamma1 C)   (M \Gamma1 C), and then taking
square roots. The lower bound is somewhat more interesting.
The matrices E 1 and E 2 are Hermitian positive semidefinite, and from the representation

it is clear that they have rank at most m and
By Corollary IV.4.9 [24], we know that
We need to show that - k are two n\Thetan matrices and the rank
of Y 2 is n\Gammam then a theorem of Weyl [20, Thm. 3.3.16] implies oe n (Y 1
Now set
I C \Gamma1
I
and notice that the eigenvalues of E 1 are the squares of the singular values of Y 1 . But
is the n \Theta n identity matrix, so by Weyl's result we obtain oe m (Y 1 ) - 1. Thus,
our conclusion follows from (12). 2
We now study the extent to which preconditioning by M mixes the signal and
noise subspaces.
Theorem 3.2. Let k be the dimension of the noise subspace and let
n\Thetam . Then
oe n\Gammak+1
Proof: Using the decompositions we have
have orthonormal columns, it follows that
oe n\Gammak+1
oe n\Gammak+1
Next we show that - oe j - oe j for oe j corresponding to the noise subspace, and thus
oe n\Gammak+1 is small. Thus, if C 1 is well-conditioned, then we are guaranteed that the
signal and noise subspaces remain unmixed.
Theorem 3.3. The (m+i)th singular value of each of the matrices C and M \Gamma1 C
lies in the interval [0; oe i
Proof: Two theorems due to Weyl for Hermitian matrices Z, Y 1 , and Y 2 with
Now from the decomposition in Equation (11), we see -n
0, and thus
Also,
We therefore likewise obtain
The proof is completed by taking square roots. 2
These theorems show that the preconditioner will be effective if C 1 is well-conditioned
and if the row sums of C
are small. We now discuss to
what extent these conditions hold for integral equation discretizations.
Property 4. Let ~
C be a Cauchy-like matrix corresponding to a real Toeplitz
matrix T that results from discretization of a smooth kernel t, normalized so that
the maximum element of T is one. Then for n sufficiently large, there exists ffl - 1
and m - n such that all elements of ~
C are less than ffl in magnitude except for those
located in four corner blocks of total dimension m \Theta m.
To understand why this is true, recall that if ~
A and ~
B are the generators of ~
where ~
A ~
, the magnitude of the (k; j)-entry of ~
C is j ~
. Thus the
largest entries in ~
appear where the numerator is large or the denominator is small.
Elements >
Elements
Elements > 3
Elements
Fig. 1. Plot revealing
The denominator of ~
C kj is j!
which is bounded above
by 2. Its smallest entries are attained for jk \Gamma jj - 0 or n, but there are very few
small values. In fact, direct computation shows that for n - 100, at least 95% of the
entries in the first row have denominators in the range [10 and the other rows
have even more in this range. Figure 1 plots values of the matrix
above given tolerance levels. As expected, there are very few large values,
and these occur only near the diagonal and the corners of the matrix.
Now consider the numerators. The formulas for A and B are determined from
direct computation in (6). The first column of A is the first unit vector, and the
second column is given by
The first column of B is
and the second column is the last unit vector. The generators for ~
C are then ~
and ~
conjugation, with F and S 0 as
described in Property 3. Therefore, the numerators are
e (n\Gamma1)
is the k th entry in the second column of ~
A and i j is the j th entry in the
first column of FS 0 B. Thus it is the normalized inverse Fourier coefficients of the
second column of A and first column of S 0 B which determine the magnitude of the
numerators, and if t is smooth, these will be large only for small indices j and k.
Therefore,
away from the corners. Thus ~
C can be permuted to contain the large elements in
the upper left block, and any pivoting strategy that produces such a permutation will
give a suitable preconditioner for our scheme.
We have observed that if Gu's algorithm is applied to a matrix with this structure,
then C 1 will contain the four corner blocks. The interested reader is referred to [10]
for details on the complete pivoting strategy, but the key fact is that Gu makes his
pivoting decisions based on the size of elements in the generator -
corresponding
to the block that remains to be factored. The resulting Cauchy-like preconditioner
C 1 for the matrix C then has the properties that the first m singular values of the
preconditioned matrix are clustered, and that the invariant subspace corresponding
to small singular values of C is not much perturbed. Thus we expect that the initial
iterations of CGLS will produce a solution that is a good approximation to the noise-free
solution.
4. Algorithmic issues. Our algorithm is as follows:
Algorithm 1: Solving
1. Compute the generators for the matrix ~
using (13)
and (14).
2. Determine an index m to define the size of the partial factorization
of ~
C and factor ~
3. Set
4. Determine the m \Theta m leading principal submatrix, C 1 ; of C and
I
5. Compute an approximate solution ~
y to M using a
few steps of CGLS.
6. The approximate solution in the original coordinate system is
y.
When to stop the CGLS iteration in order to get the best approximate solution
is a well-studied but open question (for instance, see [16] and the references therein).
We do not solve this problem, but we consider the other algorithmic issues in the
following subsections.
4.1. Determining the size of C 1 . The choice of the parameter m determines
the number of clustered singular values in the preconditioned system. It influences
the amount of work per iteration, but perhaps more importantly, the mixing of signal
and noise subspaces. We use a simple heuristic in our numerical experiments. We
compute the Fourier Transform of the data vector - g and determine the index m for
which the Fourier coefficients start to level off. This is presumed to be the noise level,
and the factorization is truncated here.
4.2. Computing the preconditioner. Since ~
C satisfies the displacement equation
(3),
are the leading principal submatrices of P
T\Omega P and Q\ThetaQ T respec-
tively, and A 1 and B 1 contain the first m rows of P T ~
A and Q T ~
respectively.
Thus the matrix C
1 has entries
~
1-i;j-n
where ~
are the elements of \Theta
and\Omega that appear in \Theta 1
respectively
and, from (5), the vectors x T
are rows of X 1 and W 1 defined as
Computing X 1 and W 1 costs O(m 2 ) operations, given the factorization of C 1 and the
matrices A 1 and B 1 .
4.3. Applying the preconditioner. Let r be a vector of length m and assume
that no pivoting was done when ~
C was factored. Heinig [17] states that C \Gamma1
may
be written as
is the jth column of X 1 , (W 1 ) j is the jth column of W 1 , and C 0 is the
Cauchy matrix C
1-i;j-m
. The notation \Delta denotes the componentwise
product of two vectors.
Fast multiplication by the matrix C 0 requires finding the coefficients of a polynomial
whose roots are the elements of \Theta 1
[6], and this process can be unstable.
To avoid this difficulty, realizing that the elements of S \Gamma1 and S 1 are roots of unity,
we extend C 0 to a matrix of size n \Theta n satisfying the displacement equation (2) with
and we develop a mathematically equivalent algorithm for
computing
Algorithm 2: Forming
For do
1. Compute
2. Extend - r by zeros so that - r is of length n.
3. Set -
r.
4. Truncate - r to length m.
5.
The product
1 r can be computed similarly.
If pivoting was done during factorization, the vector - r should be multiplied by Q
after Step 2 and by P after Step 4.
This formulation allows C
1 r to be computed in O(n lg n) operations in a stable
manner, using an observation of Finck, Heinig, and Rost [6] that any Cauchy-like
matrix can be factored as
-0.4
index
value
rhs
-0.4
index
value
solution
Fig. 2. Uncontaminated data vector (left) and exact solution vector (right) for Example 1.
minimum achieved
rel. error at iter.
43
44 2:77

Table
Minimum relative errors achieved for various values of m, Example 1.
are the Vandermonde matrices whose second columns contain
the diagonal elements
of\Omega and \Theta, respectively. The matrix H is a Hankel matrix,
i.e., one in which elements on the antidiagonals are constant. The first row is equal
to the coefficients of the polynomial
for the leading one.
Since, from Property
3,\Omega and \Theta contain roots of unity, products of the matrix C 0
with a vector are very simple to compute:
has a single non-zero diagonal extending south-west to
north-east.
the normalized, discrete, inverse Fourier Transform matrix.
is the matrix product FS 0 , where the diagonal matrix S 0 is defined in
Property 3.
Thus products C 0 - r can be computed stably in O(n lg n) operations. Since at
most, the preconditioner can be applied to a vector in O(n lg n) operations, provided
one knows X and W . This is the same order as the number of operations to apply
C to a vector, since
the product of a Toeplitz matrix with a
vector can be computed in O(n lg n) operations by embedding the matrix in a circulant
matrix [2]. Thus, each iteration of CGLS costs O(n lg n) operations.
index
value
Fig. 3. Fourier coefficients of the noisy data for Example 1.
relative
error
unpreconditioned
Fig. 4. Relative error in computed solution for Example 1.
index
value
Fig. 5. Singular values of C (solid line) and M \Gamma1 C (\Theta's) for Example 1,
5. Numerical results. In this section we summarize results of our algorithm
on three test problems using Matlab and IEEE floating point double precision arith-
metic. Our measure of success in filtering noise is the relative error, the 2-norm of
the difference between the computed estimate f and the vector -
f corresponding to
zero noise, divided by the 2-norm of -
f . In each case, we apply the CGLS iteration
with Cauchy-like preconditioner of size m. The value corresponds to no
preconditioning.
5.1. Signal processing example. As mentioned in the introduction, Toeplitz
matrices often arise in the signal processing (1-dimensional image reconstruction prob-
lems). As an example, we consider the 100 \Theta 100 Toeplitz matrix T whose entries are
defined by
ae 4
where
and
This matrix is the one used in Example 4 of [2]. The authors note that such matrices
may occur in image restoration contexts as "prototype problems" and are used to
model certain degradations in the recorded image.
The condition number of T is approximately 2:4 \Theta 10 6 . We wish to solve the
equation denotes the noisy data vector for which kek 2 =k-gk 2 , the
noise level, is about 10 \Gamma3 . The uncorrupted data, - g, and exact numerical solution, -
are displayed in Figure 2 1 . The Fourier coefficients of g are shown in Figure 3. Using
these coefficients, an appropriate cutoff value m was determined as explained in x4.1.
The solid line in Figure 4 shows the convergence of CGLS on the unpreconditioned
Toeplitz system, where the ring on the line indicates the iteration at which the minimal
value of the relative error, 2:13 \Theta 10 \Gamma1 , was achieved. Convergence of CGLS on the
preconditioned system involving the Cauchy-like matrix is also shown in Figure 4 for
two different values of m. Table 1 gives an idea of the sensitivity of the algorithm
to the choice of m, with 43 being optimal in the sense of achieving minimal
relative error among all choices of preconditioner. The number of iterations for the
preconditioned system is substantially less than for the unpreconditioned.
The singular values of T and of the preconditioned matrix M \Gamma1 C for 43 are
shown in Figure 5. As predicted by the theory in Section 3.3, the first 43 singular
values of M \Gamma1 C are clustered very tightly around one and the smallest singular values
have been left virtually untouched.
5.2. Phillips test problem. Next we consider the discretized version of the
well-known first-kind Fredholm integral equation studied by D.L. Phillips [23]. The
We first determined -
f using Matlab's square function, -
then computed -
f.
index
value
rhs
index
value
solution
Fig. 6. Uncontaminated data vector (left) and exact solution vector (right) for Example 2.
index
value
index
value
Fig. 7. Fourier coefficients of the noisy data for Example 2, two different scales.
minimum achieved
rel. error at iter.
48 4:68 \Theta 10 \Gamma2 26

Table
Minimum relative errors achieved for various values of m, Example 2.
iteration
relative
error
unpreconditioned
Fig. 8. Relative error in computed solution for 2.
index
value
Fig. 9. Singular values of C (solid line) and M \Gamma1 C (\Theta's) for Example 2,
index
value
solution
-226index
value
rhs
Fig. 10. Uncontaminated data vector (left) and exact solution vector (right) for Example 3.
kernel of the integral equation is given by t(ff; defined by
and the limits of integration are -6 and 6. We used Hansen's Matlab Regularization
Toolbox, described in [15],to generate the corresponding 400 \Theta 400 symmetric Toeplitz
matrix whose condition number was approximately 1 \Theta 10 8 . In this code, the integral
equation is discretized by the Galerkin method with orthonormal box functions. The
uncorrupted data vector is shown in Figure 6 2 . The noise level was 1 \Theta 10 \Gamma2 for this
problem.
It was difficult to determine the appropriate cutoff value m, as Figure 7 indi-
cates, but Table 2 and Figure 8 show that the savings in the number of iterations to
convergence can be substantial. In addition, for several values of m, the minimum
relative error is somewhat lower than the minimum obtained for the unpreconditioned
problem. For example, after 293 iterations, CGLS on the unpreconditioned problem
achieved a minimum relative error of 5:71 \Theta 10 \Gamma2 . For however, a minimum
relative error of 3:05 \Theta 10 \Gamma2 was reached in only 9 iterations.

Figure

9 illustrates that, as in Example 1, the first m singular values of the
preconditioned matrix are clustered around one and the singular values corresponding
to the noise subspace remain almost unchanged.
5.3. Non-symmetric example. Finally, since both previous examples involve
symmetric Toeplitz matrices, for our third example we chose to work with the 100\Theta100
matrix T whose first column is defined by
\Theta
and whose first row is h
2 We set -
f was taken to be the exact numerical solution of the problem.
Fig. 11. Fourier coefficients of the noisy data for Example 3.
minimum achieved
rel. error at iter.
28

Table
Minimum relative errors achieved for various values of m, Example 3.
iteration
relative
error
unpreconditioned
Fig. 12. Relative error in computed solution for
Fig. 13. Singular values of C (solid line) and M \Gamma1 C (\Theta's) for Example 3,
with
The condition number of T is approximately 5:31 \Theta 10 11 , making it the worst
conditioned of the three matrices. We first defined the exact solution shown in Figure
. The uncorrupted data was obtained by calculating -
f , and is also shown
in

Figure

10. White noise was added to - g to obtain the noisy data whose Fourier
coefficients are shown in Figure 11, where the noise level was determined to be 1\Theta10 \Gamma3 .
As

Figure

12 indicates, the minimum relative error obtained with no preconditioning
was 2:13 \Theta 10 \Gamma1 in 76 iterations. For values of m close to 40, however, the
preconditioned system converges in fewer than 10 iterations to the same or better
minimum relative error. We also observe from Figure 13 that in addition to clustering
the first m singular values around one, preconditioning has the benefit of reducing
the condition number.
6. Conclusions. We have developed an efficient algorithm for computing regularized
solutions to ill-posed problems with Toeplitz structure. This algorithm makes
use of an orthogonal transformation to a Cauchy-like system and iterates using the
CGLS algorithm preconditioned by a rank-m partial factorization with pivoting. By
exploiting properties of the transformation, we showed that each iteration of CGLS
costs only O(n lg n) operations for a system of n variables.
Our theory predicts that for banded Toeplitz matrices we can expect the preconditioner
determined in the course of Gu's fast modified complete pivoting algorithm
to cluster the largest singular values of the preconditioned matrix around one, keep
the smallest singular values small, and not mix the signal and noise subspaces. Thus
CGLS produces a good approximate solution within a small number of iterations.
Our results illustrate the effectiveness of our preconditioner for an optimal value of
m, and for values in a neighborhood of the optimal one. Hence, our algorithm is both
efficient and practical.
Determining the optimal value of m can be difficult, and it appears better to
underestimate the value rather than to overestimate it. Advances in computing truly
rank-revealing factorizations of Cauchy-like matrices will yield corresponding advances
in our algorithm.
Similar ideas are valid for preconditioners of the form4 C 1
are both Cauchy-like. In practice, C 2 can be determined by computing
a partial factorization of the trailing submatrix of C, remaining after C 1 is
removed. This method saves time in the precomputation of M but more iterations
may be required for convergence.
In future work, we plan to study the use of Cauchy-like preconditioners for two
dimensional problems, in which T is block Toeplitz with Toeplitz blocks, and for other
matrices related to Cauchy-like matrices.



--R

Personal Communication.
Circulant preconditioned Toeplitz least squares itera- tions
An optimal circulant preconditioner for Toeplitz systems
A lookahead Levinson algorithm for general Toeplitz systems
Regularization by truncated total least squares
An fast algorithms for Cauchy- Vandermonde matrices
Equivalence of regularization and truncated iteration in the solution of ill-posed problems
Fast Gaussian elimination with partial pivoting of matrices with displacement structure
Theory of Tikhonov Regularization for Fredholm Equations of the First Kind
Stable and efficient algorithms for structured systems of linear equations

Regularization with differential operators
Preconditioned iterative regularization for ill-posed problems
The discrete Picard condition for discrete ill-posed problems

The use of the L-curve in the regularization of discrete ill-posed problems
Inversion of generalized cauchy matrices and other classes of structured matrices

Methods of conjugate gradients for solving linear systems
Topics in Matrix Analysis
Displacement ranks of matrices and linear equations
Solution of systems of linear equations by minimized iterations
A technique for the numerical solution of certain integral equations of the first kind
Matrix Perturbation Theory
A proposal for Toeplitz matrix calculations
analysis of a fast partial pivoting method for structured matrices
The use of pivoting to improve the numerical performance of Toeplitz matrix algorithms
The rate of convergence of conjugate gradients

Pitfalls in the numerical solution of linear ill-posed problems
--TR
cauchy-like;toeplitz;conjugate gradient;least squares;preconditioner;regularization;ill-posed problems