Optimal Kronecker Product Approximation of Block Toeplitz Matrices.
This paper considers the problem of finding n  n matrices and Bk that minimize $||T - \sum A_k \otimes B_k||_F$, where $\otimes$ denotes Kronecker product and T is a banded n  n block Toeplitz matrix with banded n  n Toeplitz blocks. It is shown that the optimal and Bk are banded Toeplitz matrices, and an efficient algorithm for computing the approximation is provided.  An image restoration problem from the Hubble Space Telescope (HST) is used to illustrate the effectiveness of an approximate SVD preconditioner constructed from the Kronecker product decomposition.
Introduction
. A Toeplitz matrix is characterized by the property that its
entries are constant on each diagonal. Toeplitz and block Toeplitz matrices arise
naturally in many signal and image processing applications; see, for example, Bunch
[4] and Jain [17] and the references therein. In image restoration [21], for instance,
one needs to solve large, possibly ill-conditioned linear systems in which the coefficient
matrix is a banded block Toeplitz matrix with banded Toeplitz blocks (bttb).
Iterative algorithms, such as conjugate gradients (cg), are typically recommended
for large bttb systems. Matrix-vector multiplications can be done efficiently using
fast Fourier transforms [14]. In addition, convergence can be accelerated by preconditioning
with block circulant matrices with circulant blocks (bccb). A circulant matrix
is a Toeplitz matrix in which each column (row) can be obtained by a circular shift of
the previous column (row), and a bccb matrix is a natural extension of this structure
to two dimensions; c.f. Davis [10].
Circulant and bccb approximations are used extensively in signal and image
processing applications, both in direct methods which solve problems in the "Fourier
domain" [1, 17, 21], and as preconditioners [7]. The optimal circulant preconditioner
introduced by Chan [8] finds the closest circulant matrix in the Frobenius norm. Chan
and Olkin [9] extend this to the block case; that is, a bccb matrix C is computed to
minimize
bccb approximations work well for certain kinds of bttb matrices [7], especially
if the unknown solution is almost periodic. If this is not the case, however, the
performance of bccb preconditioners can degrade [20]. Moreover, Serra-Capizzano
and Tyrtyshnikov [6] have shown recently that it may not be possible to construct a
bccb preconditioner that results in superlinear convergence of cg.
Here we consider an alternative approach: optimal Kronecker product approxi-
mations. A Kronecker product
A\Omega B is defined as
a
Raytheon Systems Company, Dallas,
y Department of Mathematics and Computer Science, Emory University, Atlanta, GA 30322
(nagy@mathcs.emory.edu).
J. KAMM AND J. NAGY
In particular, we consider the problem of finding matrices A k , B k to minimize
s
A
where T is an n 2 \Theta n 2 banded bttb matrix, and A k , B k are n \Theta n banded Toeplitz
matrices. A general approach for constructing such an optimal approximation was
proposed by Van Loan and Pitsianis [25] (see also Pitsianis [23]). Their approach,
which we describe in more detail in Section 2, requires computing principal singular
values and vectors of an n 2 \Theta n 2 matrix related to T .
An alternative approach for computing a Kronecker product approximation T -
A\Omega B for certain deconvolution problems was proposed by Thirumalai [24]. A similar
approach for banded bttb matrices was considered by Nagy [22]. As opposed to
the method of Van Loan and Pitsianis, the schemes described in [22, 24] require
computing principal singular values and vectors of an array having dimension at most
n \Theta n, and thus can be substantially less expensive. Moreover, Kamm and Nagy [20]
show how these approximations can be used to efficiently construct approximate svd
preconditioners.
Numerical examples in [20, 22, 24] indicate that this more efficient approach can
lead to preconditioners that perform better than bccb approximations. However,
theoretical results establishing optimality of the approximations, such as in equation
(1.1), were not given. In this paper, we provide these results. In particular, we show
that some modifications to the method proposed in [22, 24] are needed to obtain an
approximation of the form (1.1). Our theoretical results lead to an efficient algorithm
for computing Kronecker product approximations of banded bttb matrices.
This paper is organized as follows. Some notation is defined, and a brief review of
the method proposed by Van Loan and Pitsianis is provided in Section 2. In Section
3 we show how to exploit the banded bttb structure to obtain an efficient scheme
for computing terms in the Kronecker product decomposition. A numerical example
from image restoration is given in Section 4.
2. Preliminaries and Notation. In this section we establish some notation to
be used throughout the paper, and describe some previous work on Kronecker product
approximations. To simplify notation, we assume T is an n \Theta n block matrix with
n \Theta n blocks.
2.1. Banded bttb Matrices. We assume that the matrix T is a block banded
Toeplitz matrix with banded Toeplitz blocks (bttb), so it can be uniquely determined
by a single column t which contains all of the non-zero values in T ; that is, some central
column. It will be useful to define an n \Theta n array P as
operator transforms matrices into vectors by stacking columns as follows:
\Theta a 1 a 2 \Delta \Delta \Delta an
a 1
a 2
an
TOPELITZ KRONECKER PRODUCT APPROXIMATION 3
Suppose further that the entry of P corresponding to the diagonal of T is known 1 .
For example, suppose that
where the diagonal of T is located at (i; is the sixth
column of T , and we write
In general, if the diagonal of T is then the upper and lower block bandwidths of
are respectively. The upper and lower bandwidths of each Toeplitz
block are
In a similar manner, the notation is used to represent a banded
point Toeplitz matrix X constructed from the vector x, where x i corresponds to the
diagonal entry. For example, if the second component of the vector
corresponds to the diagonal element of a banded Toeplitz matrix X , then
2.2. Kronecker Product Approximations. In this subsection we review the
work of Van Loan and Pitsianis. We require the following properties of Kronecker
products:
(A\Omega B)(C\Omega
(AC)\Omega (BD),
ffl If U 1 and U 2 are orthogonal matrices, then U
A more complete discussion and additional properties of Kronecker products can be
found in Horn and Johnson [16] and Graham [13].
Loan and Pitsianis [25] (see also, Pitsianis [23]) propose a general technique
for an approximation involving Kronecker products where jjT \Gamma
minimized. By defining the transformation to tilde space of a block matrix T ,
In image restoration, P is often referred to as a "point spread function", and the diagonal entry
is the location of the "point source". See Section 4 for more details.
4 J. KAMM AND J. NAGY
as
~
it is shown in [23, 25] that
s
s
(~ a k
~
where ~
a
Thus, the Kronecker product approximation
problem is reduced to a rank-s approximation problem. Given the svd of ~
it is well known [12] that the rank-s approximation ~
which minimizes jj ~
~
. Choosing ~ a
~
~ a k
~
jj F over all rank-s approximations, and thus
one can construct an approximation -
This general technique requires computing the largest s singular triplets of an
which may be expensive for large n. Thirumalai [24] and Nagy [22]
show that a Kronecker product approximation of a banded bttb matrix T can be
found by computing the largest s singular triplets of the n \Theta n array P . However, this
method does not find the Kronecker product which minimizes the Frobenius norm
approximation problem in equation (1.1). In the next section we show that if T is a
banded bttb matrix, then this optimal approximation can be computed from an svd
of a weighted version of the n \Theta n array P .
3. bttb Optimal Kronecker Product Approximation. Recall that the Van
Loan and Pitsianis approach minimizes
for a general (un-
structured) matrix T , by minimizing jj ~
(~ a k
~
k )jj F . If it is assumed that A k
are banded Toeplitz matrices, then the array P associated with the central
column of T can be weighted and used to construct an approximation which minimizes
k=1 (~ a k
~
Theorem 3.1. Let T be the n 2 \Theta n 2 banded bttb matrix constructed from P ,
is the diagonal element of T (therefore, the upper and lower block bandwidths
of T are and the upper and lower bandwidths of each Toeplitz block
are be an n \Theta n banded Toeplitz matrix with upper
lower bandwidth be an n \Theta n banded Toeplitz
matrix with upper bandwidth lower bandwidth n \Gamma j. Define a k and b k such
that A
~
~
a
~
TOPELITZ KRONECKER PRODUCT APPROXIMATION 5
s
~ a k
~
s
(W a a k )(W b b k ) T
Proof. See Section 3.1. 2
Therefore, if A k and B k are constrained to be banded Toeplitz matrices, then
can be minimized by finding a k , b k which minimize jjP
(W a a k )(W b b k ) T jj F . This is a rank-s approximation problem, involving a matrix
of relatively small dimension, which can be constructed using the svd of Pw . Noting
that W a and W b are diagonal matrices which do not need to be formed explicitly, the
construction of -
are
banded Toeplitz matrices, can be computed as follows:
ffl Define the weight vectors w a and w b based on the (i; location (in P ) of the
diagonal entry of T :
\Theta p
\Theta p
ffl Calculate
its svd
, where ": "
denotes point-wise multiplication.
ffl Calculate
a
where "./" denotes point-wise division.
The proof of Theorem 3.1 is based on observing that ~
T has at most n unique rows
and n unique columns, which consist precisely of the rows and columns of P . This
observation will become clear in the following subsection.
3.1. Proof of Theorem 3.1. To prove Theorem 3.1, we first observe that if a
matrix has one row which is a scalar multiple of another row, then a rotator can be
constructed to zero out one of these rows, i.e.,
ff
If this is extended to the case where more than two rows are repeated, then a simple
induction proof can be used to establish the following lemma.
6 J. KAMM AND J. NAGY
Lemma 3.2. Suppose an n \Theta n matrix X has k identical rows:
x
Then a sequence of k \Gamma1 orthogonal plane rotators can be constructed
such that
thereby zeroing out all the duplicate rows.
It is easily seen that this result can be applied to the columns of a matrix as well,
using the transpose of the plane rotators defined in Lemma 3.2.
Lemma 3.3. Suppose an n \Theta n matrix X contains k identical columns:
Then an orthogonal matrix Q can be constructed from a series of plane rotators such
that
\Theta p

The above results illustrate the case where the first occurrence of a row (column)
is modified to zero out the remaining occurrences. However, this is for notational
convenience only. By appropriately constructing the plane rotators, any one of the
duplicate rows (columns) may be selected for modification, and the remaining rows
(columns) zeroed out. These rotators can now be applied to the matrix ~
T .
Lemma 3.4. Let T be the n 2 \Theta n 2 banded bttb matrix constructed from P , where
ij is the diagonal entry of T . In other words,
define
~
TOPELITZ KRONECKER PRODUCT APPROXIMATION 7
Then orthogonal matrices Q 1 and Q 2 can be constructed such that
~
Proof. By definition,
representing ~
T using the n \Theta n 2 submatrices ~
~
~
~
~
it is clear that ~
contains only n unique rows, which are ~ t T
n , and that the i th
submatrix, ~
contains all the unique rows, i.e.,
~
Furthermore, it can be seen that there are occurrences of ~ t T
of ~ t T
occurrences of ~ t T
occurrences of ~ t T
of ~ t T
. Therefore, a sequence of orthogonal plane rotators can be constructed to zero
8 J. KAMM AND J. NAGY
out all rows of ~
T except those in the submatrix ~
~
W a
~
partitioning ~
~
\Theta ~
T in
where each ~
ij is an n \Theta n submatrix, it can be seen that ~
contains only n unique
columns, which are the columns of P , and that the j th submatrix ~
contains all the unique columns, i.e.,
~
Furthermore, the matrix ~
occurrences of p 1
of occurrences of pn .
Therefore, a sequence of orthogonal plane rotators can be constructed such that
~
:The following properties involving the vec and toep2 operators are needed.
Lemma 3.5. Let T , ~
T , and P be defined as in Lemma 3.4. Further, let A k be an
n \Theta n banded Toeplitz matrix with upper bandwidth lower bandwidth
and let B k be an n \Theta n banded Toeplitz matrix with upper bandwidth lower
bandwidth j. Define a k and b k such that A
Then
1. are any two matrices of the
same size,
TOPELITZ KRONECKER PRODUCT APPROXIMATION 9
2. toep2(x; are any two
vectors of the same length,
3. toep2fvec[(
4.
Proof. Properties 1 and 2 are clear from the definitions of the vec and toep2
operators. Property 3 can be seen by considering the banded Toeplitz matrices
toep(a; i) and noting that the central column of
all the non-zero entries is
an b 1
an
Therefore, property 3 holds when both sides are banded bttb matrices
constructed from the same central column, and can be extended to
applying property 2. Property 4 follows from properties 2 and 3. 2
Using these properties, Lemma 3.4 can be extended to the matrix ~
~ a k
~
k .
Lemma 3.6. Let T be the n 2 \Theta n 2 banded bttb matrix constructed from P , where
ij is the diagonal entry of T . Further, let A k be an n \Theta n banded Toeplitz matrix
with upper bandwidth lower bandwidth be an n \Theta n banded
Toeplitz matrix with upper bandwidth j. Define a k and b k such that
a
Let ~
T , W a , and W b be defined as in Lemma 3.4. Then orthogonal matrices Q 1 and
can be constructed such that
s
~ a k
~
Proof. Using Lemma 3.5,
s
A
s
a
By definition of the transformation to tilde space,
s
A
s
~
a k
J. KAMM AND J. NAGY
Applying Lemma 3.4 to T \Gamma
s
~ a k
~
a
:The proof of Theorem 3.1 follows directly from Lemma 3.6 by noting that
s
~ a k
~
s
~ a k
~
s
a
s
(W a a k )(W b b k ) T
3.2. Further Analysis. It has been shown how to minimize
when the
structure of -
T is constrained to be a sum of Kronecker products of banded Toeplitz
matrices. We now show that if T is a banded bttb matrix, then the matrix -
must adhere to this structure. Therefore, the
approximation minimizes
when T is a
banded bttb matrix.
If T is a banded bttb matrix, then the rows and columns of ~
T have a particular
structure. To represent this structure, using an approach similar to Van Loan and
Pitsianis [25], we define the constraint matrix S n;! . Given an n \Theta n banded Toeplitz
matrix T , with upper and lower bandwidths
is an
matrix such that S T
be a
4 \Theta 4 banded Toeplitz matrix with bandwidths !
TOPELITZ KRONECKER PRODUCT APPROXIMATION 11
and
Note that S T
n;! clearly has full row rank. Given the matrix T in (2.2),
~
and the rows and columns of ~
~
~
. Using the structure of ~
T , the matrix -
A
minimizing
must be structured such that A i and B i are banded Toeplitz matrices, as
the following sequence of results illustrate.
Lemma 3.7. Let
\Theta a 1 a 2 \Delta \Delta \Delta an

be the n \Theta n matrix whose structure
is constrained by S T
n;! a
be the svd of A, where
n;!
Proof. Given the svd of A,
n;!
J. KAMM AND J. NAGY
By definition, S T
n;!
Applying this result to A T , it is clear that the right singular vectors of A satisfy
if the rows of A are structured in the same manner.
Lemma 3.8. Let A =6 6 6 4
a
be the n \Theta n matrix whose structure is constrained
by S T
n;! a
i be the svd of A, where
Theorem 3.9. Let T be an n \Theta n banded block Toeplitz matrix with n \Theta n banded
Toeplitz blocks, where the upper and lower block bandwidths of T are
and the upper and lower bandwidths of each Toeplitz block are
\Theta fl u fl l
. Then
the matrices A i and B i minimizing
banded Toeplitz matrices, where the upper and lower bandwidths
of A i are given by !, and the upper and lower bandwidths of B i are given by fl.
Proof. Recall that
(~ a i
where
. The structure of T results in rank( ~
~
~
Letting ~
i be the
svd of ~
(~ a i
~
is minimized by ~
a
Therefore, A i is an n \Theta n banded Toeplitz matrix with
upper and lower bandwidths given by !, and B i is an n \Theta n banded Toeplitz matrix
with upper and lower bandwidths given by fl. 2
3.3. Remarks on Optimality. The approach outlined in this section results
in an optimal Frobenius norm Kronecker product approximation to a banded bttb
matrix. The approximation is obtained from the principal singular components of
an array Pw = W a PW b . It might be interesting to consider whether it is possible
to compute approximations which are optimal in another norm. In particular, the
method considered in [20, 22, 24] uses a Kronecker product approximation computed
from the principal singular components of P . Unfortunately we are unable to show
that this leads to an optimal norm approximation. However, there is a very close
relationship between the approaches. Since W a and W b are full rank, well-conditioned
diagonal matrices, P and Pw have the same rank. Although it is possible to establish
bounds on the singular values of products of matrices (see, for example, Horn and
Johnson [15]), we have not been able to determine a precise relationship between
the Kronecker product approximations obtained from the two methods. However we
have found through extensive numerical results that both methods give similarly good
approximations. Since numerical comparisons do not provide any additional insight
into the quality of the approximation, we omit such results. Instead, in the next
TOPELITZ KRONECKER PRODUCT APPROXIMATION 13
section we provide an example from an application that motivated this work, and
illustrate how a Kronecker product approximation might be used in practice. We
note that further comparisons with bccb approximations can be found in [20, 24].
4. An Image Restoration Example. In this section we consider an image
restoration example, and show how the Kronecker product approximations can be
used to construct an approximate svd preconditioner. Image restoration is often
modeled as a linear system:
where b is an observed blurred, noisy image, T is a large, often ill-conditioned matrix
representing the blurring phenomena, n is noise, and x is the desired true image. If
the blur is assumed to be spatially invariant, then T is a banded bttb matrix [1, 21].
In this case, the array P corresponding to a central column of T is called a point
spread function (psf).
The test data we use consists of a partial image of Jupiter taken from the Hubble
Space Telescope (hst) in 1992, before the mirrors in the Wide Field Planetary Camera
were fixed. The data was obtained via anonymous ftp from ftp.stsci.edu, in
the directory pub/stsdas/testdata/restore/data/jupiter. Figure 4.1 shows the
observed image. Also shown in Figure 4.1 is a mesh plot of the psf, P , where the
peak corresponds to the diagonal entry of T . The observed image is 256 \Theta 256, so T
is 65; 536 \Theta 65; 536.
50 100 150 200 250100200020406010305000.040.08a. Observed, blurred image. b. psf, P .
Fig. 4.1. Observed hst image and point spread function.
We mention that if T is ill-conditioned, which is often the case in image restora-
tion, then regularization is needed to suppress noise amplification in the computed
solution [21]. Although T is essentially too large to compute its condition number,
certain properties of the data indicate that T is fairly well conditioned. For instance,
we observe that the psf is not very smooth (smoother psfs typically indicate more
14 J. KAMM AND J. NAGY
ill-conditioned T ). Another indication comes from the fact that the optimal circulant
approximation of T , as well as our approximate svd of T (to be described below) are
well conditioned; specifically these approximations have condition numbers that are
approximately 20.
We also mention that if the psf can be expressed as (i.e., it has
rank 1), then the matrix T is separable. Using Theorem 3.1,
A\Omega B, where
oev). Efficient numerical methods that exploit the
Kronecker product structure of T (e.g., [2, 5, 11]) can then be used.
However, as can be seen from the plot of the singular values of P in Figure 4.2,
for this data, P is not rank one, and so T is not separable. We therefore suggest con-
values of the psf, P .
structing an approximate svd to use as a preconditioner, and solve the least squares
problem using a conjugate gradient algorithm, such as cgls; see Bj-orck [3].
This preconditioning idea was proposed in [20], and can be described as follows. Given
s
A
an svd approximation of T can be constructed as
A
. Note that the number of terms s only
affects the setup cost of calculating \Sigma. For s - 1, clearly solves the
minimization problem
min
TOPELITZ KRONECKER PRODUCT APPROXIMATION 15
over all diagonal matrices \Sigma and therefore produces an optimal svd approximation,
given a fixed
UA\Omega UB and
VA\Omega VB . This is analogous to the circulant and
bccb approximations discussed earlier, which provide an optimal eigendecomposition
given a fixed set of eigenvectors (i.e., the Fourier vectors).
In our tests, we use cgls to solve the ls problem using no preconditioner,
our approximate svd preconditioner (with terms in equation (4.1)) and the optimal
circulant preconditioner. Although we observed that T is fairly well conditioned,
we should still be cautious about noise corrupting the computed restorations. There-
fore, we use the conservative stopping tolerance jjT T

Table

4.1 shows the number of iterations needed for convergence in each case, and
in

Figure

4.3 we plot the corresponding residuals at each iteration. The computed
solutions are shown in Figure 4.4, along with the hst observed, blurred image for
comparison.

Table
Number of cgls and pcgls iterations needed for convergence.
cgls, no prec. pcgls, circulant prec. pcgls, svd prec.
43 12 4
iteration
residual
2-norm
no prec.
circulant prec.
svd prec
Fig. 4.3. Plot of the residuals at each iteration.
5. Concluding Remarks. Because the image and psf used in the previous
section come from actual hst data, we cannot get an analytical measure on the accuracy
of the computed solutions. However, we observe from Figure 4.4 that all
solutions appear to be equally good restorations of the image, and from Figure 4.3
we see that the approximate svd preconditioner is effective at reducing the number
of iterations needed to obtain the solutions. Additional numerical examples comparing
the accuracy of computed solutions, as well as computational cost of bccb and
the approximation svd preconditioner, can be found in [19, 20]. A comparison of
J. KAMM AND J. NAGY
50 100 150 200 25010020050 100 150 200 250100200a. hst blurred image. b. cgls solution, 43 iterations.
50 100 150 200 25010020050 100 150 200 250100200c. pcgls solution, circ. prec., 12 its. d. pcgls solution, svd prec., 4 its.
Fig. 4.4. The observed image, along with computed solutions from cgls and pcgls.
computational complexity between bccb preconditioners and the approximate svd
preconditioner depends on many factors. For example:
ffl What is the dimension of P (i.e., the bandwidths of T )?
ffl Is a Lanczos scheme used to compute svds of P , A 1 and
ffl Do we take advantage of band and Toeplitz structure when forming matrix-matrix
products involving UA , UB , VA , VB and A k , B k ,
TOPELITZ KRONECKER PRODUCT APPROXIMATION 17
ffl How many terms, s, do we take in the Kronecker product approximation?
ffl For bccb preconditioners: is n a power of 2?
If we assume T is set up and application of the approximate
svd preconditioner is at most O(n 3 ). If we further assume that n is a power of 2, then
the corresponding cost for bccb preconditioners is at least O(n 2 log 2 n). It should be
noted that the approximate svd preconditioner does not require complex arithmetic,
does not require n to be a power of 2, or any zero padding. Moreover, decomposing
T into a sum of Kronecker products, whose terms are banded Toeplitz matrices,
might lead to other fast algorithms (as has occurred over many years of studying
displacement structure [18]). In this case, the work presented in this paper provides
an algorithm for efficiently computing an optimal Kronecker product approximation.



--R


Restoration of images degraded by spatially varying pointspread functions by a conjugate gradient method

Stability of methods for solving Toeplitz systems of equations
Application of ADI iterative methods to the image restoration of noisy images
Any circulant-like preconditioner for multilevel Toeplitz matrices is not superlinear
Conjugate gradient methods for Toeplitz systems
An optimal circulant preconditioner for Toeplitz systems
Preconditioners for Toeplitz-block matrices

Algorithms for the regularization of ill-conditioned least squares problems with tensor product structure
Matrix Computations
Kronecker Products and Matrix Calculus: with Applications
Restoration of atmospherically blurred images by symmetric indefinite conjugate gradient techniques
Matrix Analysis

Fundamentals of Digital Image Processing
Theory and applications
Singular value decomposition-based methods for signal and image restoration
Kronecker product and SVD approximations in image restoration
Iterative Identification and Restoration of Images
Decomposition of block Toeplitz matrices into a sum of Kronecker products with applications in image restoration
The Kronecker Product in Approximation and Fast Transform Generation
High performance algorithms to solve Toeplitz and block Toeplitz matrices
Approximation with Kronecker products
--TR

--CTR
S. Serra Capizzano , E. Tyrtyshnikov, How to prove that a preconditioner cannot be superlinear, Mathematics of Computation, v.72 n.243, p.1305-1316, July
conjugate gradient method;block Toeplitz matrix;singular value decomposition;kronecker product;image restoration;preconditioning