Numerical Condition of Discrete Wavelet Transforms.
The recursive algorithm of a (fast) discrete wavelet transform, as well as its generalizations, can be described as repeated applications of block-Toeplitz operators or, in the case of periodized wavelets, multiplications by block circulant matrices. Singular values of a block circulant matrix are the singular values of some matrix trigonometric series evaluated at certain points. The norm of a block-Toeplitz operator is then the essential supremum of the largest singular value curve of this series. For all reasonable wavelets, the condition number of a block-Toeplitz operator thus is the lowest upper bound for the condition of corresponding block circulant matrices of all possible sizes. In the last section, these results are used to study conditioning of biorthogonal wavelets based on B-splines.
Introduction
. Orthogonality is a very strong property. It might exclude
other useful properties like, for example, symmetry in the case of compactly supported
wavelets [6, 7]. Consequently, in many applications biorthogonal wavelets have
been used rather than the orthogonal ones. Stability of such bases has been studied
and conditions for Riesz bounds to be finite were established [2, 3, 4, 5]. However,
when dealing with applications, one would like to have some quantitative information
about sensitivity to such things like noise in the data or quantization. Some relevant
estimates can be found in the engineering literature on multirate filter banks, where
noise is modelled as a random process and its transmission through the system is
studied; see, e.g., [12]. However, most of these results concern particular designs and
implementations. Here we will use an alternative approach-we will look at discrete
wavelet transforms from the point of view of linear algebra.
For example, let us consider the process of image compression using wavelets
(see, e.g., [1, 11, 14]). The algorithm has three steps. First, the discrete wavelet
transform is applied to the image, then the resulting data is quantized and finally
it is coded in some efficient way. The purpose of the transform is to increase the
compressibility of the data and to restructure the data so as, after decompression,
the error caused by quantizing is less disturbing for a human viewer then if the image
was quantized directly without a transform. The encoded image can be manipulated
in different ways (e.g., transmitted over network) which can cause further distortions.
To decompress the image we just need to decode the data and to apply the inverse
transform. Let us denote the error vector that is added to the transformed data y
before the reconstruction by u and let us suppose that we know the magnitude of the
denotes the original image, the relative error
in the reconstructed image is
This work was supported by the Flinders University of South Australia, the Cooperative Research
Centre for Sensor Signal and Information Processing, Adelaide, and the Australian Government
y National Institute of Standards and Technology, Gaithersburg, MD 20899-0001
A
If no further assumptions are imposed on the image and type of the error, this estimate
is the best possible. Also in other applications, the sensitivity to errors can be shown
to be naturally related to the condition number of the transform matrix with respect
to solving a system of linear equations,
Condition number depends on the norm. For finite matrices we will use here
matrix 2-norm, which is induced by the Euclidean vector norm. When necessary, we
will use subscript 2 to emphasise that we deal with these norms. We will speak also
about condition number of an operator l 2 ( Z
Z ). We define it also by (1.1);
the norm is the operator norm induced by the norm of l 2 ( Z
Due to the translational character of wavelet bases, matrices and operators involved
happen to have a characteristic structure-they are block circulant and block-
Toeplitz, respectively. This structure can be employed when the condition numbers
are computed; Fourier techniques can be used to transform them to a block diagonal
form. This then leads to studying the (point-wise) singular values of certain trigonometric
matrix series. In (Section 3 we study the finite case. The singular values of
a block circulant matrix are shown to be the singular values of small matrices arising
from the "block discrete Fourier transform" of the first block row of the block
circulant matrix. In (Section 4 we generalize this result for block-Toeplitz operators
Z ). Situation is rather more complicated there, because the Fourier
transform maps the discrete space l 2 ( Z
Z ) onto the functional space L 2 ([0; 2-)). As the
main result we show there that
ess sup
oe
where C(A) is the block-Toeplitz operator the infinite matrix of which is generated
by strip
Z , being square blocks) and
The proof is based on the point-wise singular decomposition of A; some difficulties
arising from the fact that we have to ensure that the singular vector we want to construct
has square integrable components must be overcome on the way. For reasonable
wavelets, the curves of singular values of A have some smoothness and essential supremum
and infimum become supremum and infimum or even maximum and minimum.
Condition number of C(A) is then lowest upper bound on condition of periodized
wavelet transforms for all possible lengths of data. We also describe how some particular
properties of the wavelets imply a certain structure of the singular values. These
observations can be used to further improve the efficiency of computing the condition
numbers.
In the last section of this paper, we apply this technique to study conditioning
of biorthogonal B-spline wavelets constructed by Cohen, Daubechies and Feauveau
[5], nowadays probably the most often applied biorthogonal wavelets. We show there
NUMERICAL CONDITION OF DISCRETE WAVELET TRANSFORMS 3
that the condition number increases exponentially with the order of the spline. Conditioning
can be significantly improved by suitable scaling of the wavelet functions,
but, even for the optimal scaling, the growth has exponential character.
After finishing the first version of this paper, I became familiar with related
works by Keinert [8] and Strang [9]. While Strang's work concerns mostly Riesz
bounds for subspaces in a multiresolution analysis and wavelet decomposition, Keinert
concentrates on conditioning of finitely sized transforms and asymptotic estimates for
deep recursive transforms. He also presents a number of numerical experiments that
show how these estimates are realistic when some specific types of introduced errors
are considered (e.g., white noise). In this revised version I have tried to emphasise
results that are complementary to those of Keinert and Strang.
2. Translational and wavelet bases and the operators of the change
of a basis. Let us consider some translation-invariant subspace of L 2 ( R
I ) with a
translational Riesz basis fu k
Zg generated by some r-tuple
of functions u being the translation step. Let this subspace have
another, similar, basis fv k
Z g. Each of the functions v k ,
can be expressed in the terms of the first basis; there exist sequences
Z ) such that
r
a (k;l)
Let us form from these coefficients r \Theta r matrices An , n 2 Z
n will be the element
of An in kth row and lth column. We denote by A the infinite strip of concatenated
matrices An , n 2 Z
Z ,
and we define C(A) to be an infinite block-Toeplitz matrix
We will denote by C(A) also an operator l 2 ( Z
Z ) that can be represented by
such a matrix. If
r
ff
r
for some l 2 ( Z
to
fff n gn2Z Z , that is, it is the operator of the change of a basis.
Because of practical reasons (handling of finite data), periodized bases are often
used in the wavelet context. If, for some integer N , we denote
per
A
then fu per
are bases for some subspace of L 2 ([0; Nh)) and the operator of the
change of basis from the latter to the former can be represented by a block circulant
matrix
SN
where
A multiresolution analysis is a sequence of embedded subspaces of L 2 ( R
I ) generated
by the translates of an appropriately dilated scaling function. In particular,
Z g:
There are wavelet subspaces generated by a wavelet function,
Z
and these subspaces satisfy
The scaling and wavelet function thus have to conform to the two-scale relations that
are usually written as
In the (fast) discrete wavelet transform, we perform recursively the change of basis
from f2 j=2 '(2
Z g to f2 (j \Gamma1)=2
. We can consider both bases to be generated
by two functions, former by u 1), the latter
by x). The translation step h is
here. The recursive inverse transform thus can be associated with repeated
applications of C(A)   , where
that is,
In fact, C (
A and the recursive transform itself can be seen as
repetitive applications of a block-Toeplitz operator. As in the case of A,
~
NUMERICAL CONDITION OF DISCRETE WAVELET TRANSFORMS 5
sequences f ~ hn gn2Z Z and f~g n gn2Z Z determine the biorthogonal counterparts of the
scaling and wavelet function, ~
' and ~
/ by relations analogous to (2.1).
Although the conditioning of this basic step of recursive transform is crucial, we
want to study also how the error cumulates in the recursive transform. Since all the
bases involved have translational character, we can use the same approach as for one
step for the transform of any finite depth; we can always find a common translation
step. For example, let us consider two steps of recursion. We perform, in fact the
change of basis in V j from f2 j=2 '(2
Z g to f2 (j \Gamma2)=2
Z
Z g. All these
bases can be considered to be translational bases with translation step
generated by four functions. We have
and
The infinite strip A thus will have four rows; the entries can be easily found by
recursive applications of (2.1). In particular, if we denote the sequences that form
rows of A by fb
0 , we have
b (1)
b (2)
An analogous approach can be used for generalizations of classical wavelet transforms
like those based on more than one scaling and wavelet function and general
integer dilation parameter m - 2 (multiwavelets, higher multiplicity wavelets) or
non-stationary wavelets, where different block-Toeplitz operators applied in the recursive
algorithm. Also wavelet packets transforms, where also wavelet spaces are
further decomposed, can be described in a similar way.
3. Numerical condition of block circulant matrices. Any circulant matrix
is unitarily similar to a diagonal matrix. This matrix has (up to scale) the discrete
Fourier transform of the first row of the original matrix on the diagonal and the
similarity matrix is the matrix of the discrete Fourier transform itself. This fact can
be generalized for block circulant matrices as follows.
Theorem 3.1. Each block circulant matrix is unitarily similar to a block diagonal
matrix. In particular, CN (A) is similar to a matrix with diagonal blocks equal to
Proof. Let !N is the primitive Nth root of unity, us first
create the matrix of the "block discrete Fourier transform";
I
I
I
(3.
6 R. TURCAJOV '
A
I being the r \Theta r identity matrix. Such a matrix is unitary and the r \Theta r block in
row and (n + 1)th block column
of\Omega r;N CN
(A)\Omega
An\Gammal+kN
being the Kronecker delta.
Since the singular values are preserved by unitary transformations and the singular
values of a block diagonal matrix are the singular values of the diagonal blocks,
the theorem above has the following corollary.
Corollary 3.2. A number oe is a singular value of CN (A) if and only if it is a
singular value of A(2-in=N) for some
Let us remind here, that the 2-norm of a matrix M equals to its largest singular
value, which we will denote here oe max (M ). Similarly, oe min (M) will stand for the
smallest singular value, the 2-norm of M \Gamma1 .
Corollary 3.3.
If N 1 is a divisor of N , cond 2 (CN1 (A)) - cond 2 (CN (A)), because all the singular
values of CN1 (A) are simultaneously singular values of CN (A). This means that, for
the recursive transform, we could estimate the condition in each step by the condition
number of the largest block circulant matrix involved, applied in the first step of the
recursion, since in each next step just m-times smaller matrix is used, m being the
dilation factor.
It would be useful to have some estimate completely independent of the size of
the block circulant matrix. One such estimate is straightforward,
Notice that if the curves of the largest and smallest singular values are continuous
(which happens, for example, for compactly supported wavelets, when A contains only
a finite number of non-zero entries) this is the lowest upper bound for cond(CN (A))
independent of N . We will show in the next section that for any reasonable wavelet
the right hand side of (3.2) represents, in fact, the condition number of C(A).
4. Norm and condition number of block-Toeplitz operators. Similarly
as in the previous section, we will apply here a "block Fourier transform". However,
here the situation is a little more complicated than in the case of finite matrices.
Let us denote l 2
Z ) the Hilbert space of (column) vectors of length r with all
components in l 2 ( Z
Z ). We can see this space also as a space of vector-valued sequences.
The inner product is
r
r
a
NUMERICAL CONDITION OF DISCRETE WAVELET TRANSFORMS 7
subscripts determine entries of sequences, while superscripts entries of vectors. Sim-
ilarly,
r ([0; 2-)) is the Hilbert space of r-vectors of square integrable functions on
[0; 2-) with the inner product
r
r
Z 2-f
Z 2-g(-)   f(-) d-:
To find the norm of the operator C(A) induced by the norm of l 2 ( Z
Z ), we employ
Hilbert space isomorphisms of these spaces. First, there is a trivial isomorphism
between l 2 ( Z
Z ) and l 2
Z
Second, component-wise Fourier transform is a Hilbert space isomorphism l 2
r ([0; 2-)). For a sequence c 2 l 2 ( Z
Z ) the Fourier transform bc 2 L 2 ([0; 2-)) is defined
as
where the sum converges in L 2 ([0; 2-)) sense. Since 1
e \Gammaik- , k 2 ZZ is an orthonormal
basis for L 2 ([0; 2-)), the inverse mapping is given by
Z 2-bc(-)e ik- dand
the Fourier transform as defined above is a Hilbert space isomorphism l 2 ( Z
2-)). The extension to the vector case is obvious.
Infinite Toeplitz matrices represent convolution operators. For sequences a; b 2
Z ), the convolution c = a   b has entries
Convolution operators are closely related to multipliers. The link is the Fourier transform

Lemma 4.1. Let a; b 2 l 2 ( Z
Z ) and let a   b 2 l 2 ( Z
Z ) or ba b b 2 L 2 ([0; 2-)). Then
d
a
Proof. For any l 2 Z
Z ,
d
Because the Fourier transform is a Hilbert space isomorphism,
A
The last term represents the lth entry of the inverse Fourier transform of
Theorem 4.2. The operator l 2 ( Z
represented by C(A) is isomorphic
with a matrix multiplier
k2ZZ A k e ik- that maps L 2
r ([0; 2-)) \Gamma! L 2
r ([0; 2-)),
u(-) \Gamma! A(-)u(-).
Proof. By the former isomorphism, C(A) is isomorphic with the operator l 2
l 2
Z ), for which d, the image of c, is given by the formula
A k\Gammal c k ; l 2 Z
We will slightly abuse the notation and denote this operator also by C(A).
Since we assume that C(A) represents the change from one Riesz basis to another,
Z ) and the series
k2ZZ A k e ik- converges component-wise in
straightforward calculation shows that (4.1) can be extended to
matrix/vector case (the Fourier transform being defined component-wise). Because
fA \Gammak g k2ZZ
A(\Gamma-),
d
A(\Gamma-)b c(-);
Z ) or b
r ([0; 2-)). A convolution-type operator
thus becomes in the Fourier domain, indeed, the matrix multiplier A.
The norm of C(A) induced by l 2 ( Z
thus equals to the norm of the matrix multiplier
A as an operator L 2
r ([0; 2-)) \Gamma! L 2
r ([0; 2-)). The following theorem gives
formulae for the norm of a multiplier and its inverse.
Theorem 4.3. Let A be an r \Theta r matrix multiplier with measurable components.
Then
sup
kAuk ess sup
r ([0;2-))
kAuk ess inf
min
Proof. Let us set
ess sup
and let x 2 L 2
r ([0; 2-)). Then
r ([0;2-))
and hence
sup
r ([0;2-))
r ([0;2-
We need to show that we have the lowest upper bound, in other words, that for each
NUMERICAL CONDITION OF DISCRETE WAVELET TRANSFORMS 9
Let us take point-wise the singular value decomposition of A;
where V(-) and U(-) are r \Theta r unitary matrices and \Sigma(-) is the diagonal matrix with
the singular values of A(-) on the diagonal, in decreasing order. We will denote these
singular values oe j (-), To construct x ffl we need a path of right singular
vectors corresponding to the largest singular value, something like the first column of
U , but we have to ensure that this path is square integrable.
First, oe 1 is a measurable function, because A has measurable components
and the matrix norm is a continuous function of the entries. Let us define
Then C has measurable components and, for k ! +1, C(-) k \Gamma! P(-), where
and D(-) is a diagonal matrix with the elements on the diagonal equal to either 1 or 0;
if oe 1 (-) is of multiplicity m (m depending on -), then first m elements are 1 and all the
others are 0. Notice that P(-) is the orthogonal projector onto the subspaces spanned
by all right singular vectors corresponding to singular values oe 1
Because P is the limit of a sequence of matrices with measurable components, its
components are measurable too.
Now, for any ffl ? 0, the set
is a measurable set and -(S ffl ) ? 0. Since P(-) 6= 0 for any -, there exist j and a set
~
that p(-), the jth column of P(-), is non-zero for - 2 ~
Let us set
-( ~
Because x ffl has measurable components and j~x
have
r ([0; 2-)). A simple calculation shows that We have
consequently,
and
Finally,
-( ~
A
which finishes the proof of the first part of the theorem.
us concentrate on the second part of the statement. Let us denote
~
ess inf
min
Clearly, for any u, kuk
kAuk
Z 2-'
min
We now have to show that for every ffl ? 0 there exists x ffl ,
that
In order to do that we first need to construct a square integrable path of right singular
vectors corresponding to the path of the smallest singular values, oe r .
Let us take, again, point-wise singular value decomposition of A,
Now, for a positive integer k, let us consider a matrix A(-)   A(-)
k I . We have
I
therefore such a matrix is invertible and the norm of the inverse is (oe r
If we set
oe r (-) 2 +k
I
then C k has measurable components and C k (-) l \Gamma! P(-), k ! +1, l ! +1, where
is, again, a diagonal matrix with the elements on the diagonal equal to either 1
or 0, but now, if oe r (-) is of multiplicity m, then the first r \Gamma m elements are 0 and
all the others are 1. The components of the matrix P are measurable functions and
the matrix P(-) is now the orthogonal projector onto the subspace spanned by right
singular vectors corresponding to the singular values oe r\Gammam+1
for any vector x(-) of unit norm from its range,
The rest of the proof would follow the lines of the proof of the first part, with S ffl being
chosen as
NUMERICAL CONDITION OF DISCRETE WAVELET TRANSFORMS 11
we would have
oe r (-) 2
-( ~
The norm of A induced by the norm of L 2
r ([0; 2-)) thus is
ess sup
oe
the mapping is invertible if and only if
ess inf
min
and the norm of the inverse equals
ess inf
oe min (A(-)):
Combining the results above we obtain the following theorem.
Theorem 4.4. The condition of the operator C(A) (in the norm induced by the
norm of l 2 ( Z
Z
ess sup -2[0;2-) oe max (A(-))
ess inf -2[0;2-) oe min
For all wavelets of practical interest, A has only a finite number of non-zero
entries or at least the sequences forming its rows decay very fast. This implies some
smoothness of entries of A and, consequently, essential supremum of oe max and essential
infimum of oe min coincide with the supremum and infimum, respectively. As we already
pointed out, cond(C(A)) then represents sup N cond(CN (A)).
Let us make a few comments about the structure of singular values of A(-) in
relation with some special properties of A. First, when the underlying bases comprise
of real functions, the entries of A are real and, consequently,
means that the singular values in - \Gamma - and - are the same and we can restrict
our attention onto interval [0; -], only.
Another interesting effect is caused by all the scaling and wavelet functions and
their biorthogonal counterparts being compactly supported. This corresponds to the
fact that only a finite number of square blocks both in A and in ~
A that generates
are non-zero. It is well known, particularly in the filter bank
context (see, e. g., [12], [13]), that this happens if and only if there exist a non-zero
constant ff and an integer p such that
det(
A k z \Gammak
for any z 2 CI , z 6= 0. Because determinant is the product of singular values, the
equation above implies that
Y
A
for some positive constant fi independent of -. This is particularly useful when A has
only two rows. The singular values of A(-) are then inversely proportional and the
maximum and minimum over - then occur at the same point. That is,
5. Alternative expression. Let the sequences that form the rows of A be
. Sometimes it is easier to deal with Fourier series
than with A. We will see an example in (Section 6, when we will study conditioning
of biorthogonal spline wavelets. In these cases it is better to use a different matrix
function.
Theorem 5.1. A number oe is a singular value of A(-) if and only if it is a
singular value of B(\Gamma-=r), where
b (2) (-) b (2) (-
.
Proof. Using the notation introduced in the proof of (Theorem 3.1, we have, for
any
2-l
r
r =r
r
a (s;k+1)
This is because
r equals to r if divisible by r and it is 0 otherwise.
Consequently,
B(-)\Omega
rA(\Gammar-)D r (-);
where\Omega 1;r is the r \Theta r matrix of the discrete Fourier transform and D r (-) is the
diagonal matrix with the diagonal entries equal to e \Gammaik- ,
particular order).
1;r is unitary and so is D r (-) (for any -), the statement of
the theorem holds.
Just let us point here that, instead of considering each row of A separately, we
could use block rows, each of them comprising of, let say, p rows. We then would
obtain similar result with some p \Theta p matrices B
n instead of scalars b
instead of
\Omega 1;r we would
use\Omega p;r=p and, similarly, D r (-) would be replaced by a matrix with
diagonal blocks equal to e \Gammaik- I , This might be useful for the
case of multiwavelets (more than one scaling function) when the two-scale equations
analogous to (2.1) have matrix coefficients (see, e.g., [10]).
NUMERICAL CONDITION OF DISCRETE WAVELET TRANSFORMS 13
6. Conditioning of biorthogonal wavelets based on B-splines. B-splines.
Biorthogonal wavelets based on B-splines were introduced by Cohen, Daubechies and
Feauveau in [5]. To the B-spline basis function of a particular order (which represents
a scaling function) there exists a whole family of possible biorthogonal counterparts
with different size of support and regularity. We will use here the notation of [2],
where the sequences determining the scaling and wavelet functions through the two-scale
equations of type (2.1) are given in terms of trigonometric polynomials
~
Since the scaling function ' equals to the B-spline of order n,
For any integer K such that 2K -
~
determines a biorthogonal scaling function; PK is the solution of Bezout problem
in particular,
The
corresponding to the wavelet filters are
then defined as
We have
and, because we deal with compactly supported real classical wavelets with dilations
by 2, we are interested in the maximum of the condition number of B(-) on [0; -=2].
Squares of the singular values of B(-) are the eigenvalues of the matrix B(-)B(-)
and they satisfy a quadratic equation
Fairly straightforward although somewhat tedious calculations show that the coefficients
of this equation are
det(B(-)B(-)
14 R. TURCAJOV '
A
and
where
Theorem 6.1. The numerical condition of one level of the (fast) discrete wavelet
transform based on B-spline biorthogonal wavelets of order n defined above is at least
independently of the value of K.
Proof. Since
(cf. (6.1)), substituting -=2 into the formulae above we obtain
Squares of the singular values of B(-=2) thus equal \Gamman , respectively, and the
condition of this matrix is 2 n . The condition hence must be at least 2 n .
Numerical experiments show that the condition number often equals 2 n . From
the point of view of conditioning, it is better to choose K smaller for low order splines
and larger for higher order splines; see Table 1 at the end of the paper.
Once the scaling filters m 0 and ~
are given, (6.2) is not the only possibility for
the corresponding wavelet filters. The entire freedom can be described as follows:
Z , ff 6= 0. The choice of k is, from the point of view of numerical condition,
irrelevant, but the scaling by ff can be used to improve the condition. In the case
of the spline wavelets improvement can be significant. However, it turns out that
whatever scaling we choose, we can't beat the exponential growth with the order of
the spline.
Theorem 6.2. For any scaling factor ff, the condition of one step of discrete
wavelet transform with spline biorthogonal wavelet of order n is at least 2 n=2 .
Proof. Instead of the condition of B(-) we need to study here the condition of
before, in (6.2). For -=2
the singular values of B ff (-=2) are jffj 2 n=2 and 2 \Gamman=2 and its condition hence is jffj 2 n
for \Gamman and 1=(jffj2 n ) for \Gamman . On the other hand, for
and its condition is jffj for jffj - 1 and 1=jffj for jffj ! 1. Combining these results
we see that the condition of the wavelet transform can not be better than jffj2 n if
NUMERICAL CONDITION OF DISCRETE WAVELET TRANSFORMS 15
Consequently, whatever jffj we choose, the
condition is at least 2 n=2 .
The optimal scaling parameter is usually equal or close to 2 \Gamman=2 , see Table 2 and
3. Notice that this is true especially for the wavelets that have condition number
equal to 2 n . The condition of the optimally scaled wavelet then equals 2 n=2 , in most
cases.

Figures

1-5 show some typical behaviour of the singular value curves in dependence
on the order of the spline, parameter K, scaling of the wavelet and depth of the
transform. There are some interesting details there like, for example, the presence of
points where the plot looks almost like if two curves were intersecting each other, but,
in fact, we have two different curves that have turning points and are well separated.

Acknowledgements

. The author was a postgraduate research scholar supported
by the Australian government. She thanks Jaroslav Kautsky for suggesting the topic
and many fruitful discussions.



--R

Image coding using wavelet transform
in Wavelets: A Tutorial in Theory and Applications
A stability criterion for biorthogonal wavelet bases and their related subband coding scheme

Biorthogonal bases of compactly supported wavelets
Orthonormal bases of compactly supported wavelets

Numerical stability of biorthogonal wavelet transforms
Inner products and condition numbers for wavelets and filter banks
Short wavelets and matrix dilation equations
Compact image coding using wavelets and wavelet packets based on non-stationary and inhomogeneous multiresolution analyses
Multirate Systems and Filter Banks
Perfect reconstruction FIR filter banks: Some properties and factorizations

--TR

--CTR
Zhong-Yun Liu, Some properties of centrosymmetric matrices, Applied Mathematics and Computation, v.141 n.2-3, p.297-306, 5 September
block-Toeplitz operators;translational bases;biorthogonal wavelets;numerical condition;block circulant matrices