Sparse Image Representation Techniques
Sparse Image Representation Techniques
a dissertation
submitted to the department of statistics
and the committee on graduate studies
of stanford university
in partial fulfillment of the requirements
for the degree of
doctor of philosophy
Xiaoming Huo
August 1999
c Copyright by Xiaoming Huo 1999
All Rights Reserved
ii
I certify that I have read this dissertation and that in my
opinion it is fully adequate, in scope and quality, as a disser-
tation for the degree of Doctor of Philosophy.
David L. Donoho
(Principal Adviser)
Stephen Boyd
Michael A. Saunders
iii
iv
To my mother and father
and my sister, Xiao-bin
v
vi
To find a sparse image representation.
vii
viii
Abstract
ix
x
Acknowledgments
I would like to thank my advisor, Professor David L. Donoho, for introducing me to the
areas of statistical signal/image processing and optimization, and for providing me with a
continuous stream of suggestions, feedbacks and encouragements during the period I worked
on this thesis. He is the best scientist I have ever seen, and his penetrating views on many
scientific fields have kept surprising me.
I am indebted to Professor Michael Saunders for letting me use and modify his soft-
ware, and for helping me in programming and technical writing. His expertise in scientific
computing is something I have always admired.
I would like to thank Professor Stephen Boyd for giving a tremendous amount of sug-
gestions during the group meetings and for serving on my reading committee. His clarity
and style in technical presentation and sharpness in abstract thinking are something that I
have always tried to imitate.
I would like to thank Professor Brad Efron and Professor Iain Johnstone for writing
recommendation letters on my behalf during my job search. I would also like to thank
Professor Jim Dai for helping me create the interview opportunity in Georgia Institute of
Technology, which will soon be the next stop in my professional career.
I want to thank all the members of the Department of Statistics at Stanford for creating
a stimulating and supporting environment. I would like to thank all the members in the
Boyd-Donoho-Saunders group, whose presence sometimes is my major source of momentum.
I want to thank all my friends at Stanford, who have made another part of my life full of
joy and excitement. (I am afraid to list all their names here as we usually do, because I
know that no matter how hard I try, the list will always be incomplete.)
Finally, I want to express my deepest gratitude to my parents and my sister for constant
and unconditional love and support. Without them, none of my accomplishments would be
possible. To them, I dedicate this thesis.
xi
xii
Contents
Abstract ix
Acknowledgments xi
Content xvii
Nomenclature xxv
1 Introduction 1
1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
xiii
2.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4 Proof . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
xiv
4.7 Newton Direction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.8 Comparison with Existing Algorithms . . . . . . . . . . . . . . . . . . . . . 93
4.9 Iterative Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
4.10 Numerical Issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.11 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.11.1 Connection With Statistics . . . . . . . . . . . . . . . . . . . . . . . 96
4.11.2 Non-convex Sparsity Measure . . . . . . . . . . . . . . . . . . . . . . 97
4.11.3 Iterative Algorithm for Non-convex Optimization Problems . . . . . 97
4.12 Proofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.12.1 Proof of Proposition 4.1 . . . . . . . . . . . . . . . . . . . . . . . . . 98
4.12.2 Proof of Theorem 4.1 . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4.12.3 Proof of Theorem 4.2 . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.12.4 Proof of Theorem 4.3 . . . . . . . . . . . . . . . . . . . . . . . . . . 102
xv
6 Simulations 119
6.1 Dictionary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
6.2 Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
6.3 Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.4 Decay of Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
6.5 Comparison with Matching Pursuit . . . . . . . . . . . . . . . . . . . . . . . 125
6.6 Summary of Computational Experiments . . . . . . . . . . . . . . . . . . . 127
6.7 Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
6.8 Scale of Efforts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
xvi
B.2.3 X-interpolation: from Cartesian to Polar Coordinate . . . . . . . . . 154
B.2.4 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
B.3 Adjoint of the Fast Transform . . . . . . . . . . . . . . . . . . . . . . . . . . 159
B.4 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
B.4.1 Storage and Computational Complexity . . . . . . . . . . . . . . . . 161
B.4.2 Effective Region . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
B.4.3 Ill-conditioning of the Fast X-interpolation Transform . . . . . . . . 162
B.5 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
B.5.1 Basic Elements for the Fast Edgelet-like Transform . . . . . . . . . . 163
B.5.2 Edgelet-like Transform for Some Artificial Images . . . . . . . . . . . 164
B.5.3 Edgelet-like Transform for Some Real Images . . . . . . . . . . . . . 165
B.6 Miscellaneous . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
B.6.1 Examples of Interpolation Functions . . . . . . . . . . . . . . . . . . 167
Bibliography 170
xvii
xviii
List of Tables
3.1 Comparison of transforms and the image features that they are good at
processing. The third column lists the order of computational complexity for
their discrete fast algorithms. . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2 Definitions of four types of DCTs. . . . . . . . . . . . . . . . . . . . . . . . 39
3.3 Number of multiplications and additions for various fast one-D DCT/IDCT
algorithms. The number in { · } is the value when N is equal to 8. . . . . 42
3.4 Decomposition of covariance matrix that can be diagonalized by different
types of DCT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
xix
xx
List of Figures
2.4 Diagram of image transmission using transform coding. The sequence (or
vector) y is the coefficient sequence (or vector) and ỹ is the recovered coeffi-
cient sequence (or vector). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.5 Intuition to compare the speeds of decay for two sequences. A sequence whose
sorted-amplitudes-curve B is, as we consider, sparser than the sequence whose
sorted-amplitudes-curve is A. . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1 Sampling of four types of DCT at the index domain. For example, both the
first index and the second index for the type-I DCT take integral values. . . 40
xxi
3.4 Illustration of the discrete algorithm for forward orthonormal wavelet trans-
form on a finite-length discrete signal. The upper graph is for cases having
3 layers. The width of each block is proportional to the length of the corre-
sponding subsequence in the discrete signal. The bottom one is a symbolic
version for general cases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
3.5 Multiresolution analysis of a point singularity with Haar wavelets. . . . . . 62
3.6 Two-dimensional wavelet basis functions. These are 32 × 32 images. The
upper left one is a tensor product of two scaling functions. The bottom right
2 by 2 images and the (2, 2)th image are tensor products of two wavelets.
The remaining images are tensor products of a scaling function with a wavelet. 63
3.7 Idealized tiling on the time-frequency plane for (a) sampling in time domain
(Shannon), (b) Fourier transform, (c) Gabor analysis, (d) orthogonal wavelet
transform, (e) chirplet, (f) orthonormal fan basis, (g) cosine packets, and (h)
wavelet packets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
A.1 Edgelet transform of the Chinese character “Huo”: (a) is the original; (d) is
the sorted coefficients; (b), (c), (e) and (f) are reconstructions based on the
largest 200, 400, 800, 1600 coefficients, respectively. . . . . . . . . . . . . . . 138
A.2 Edgelet transform of the sticky image: (a) is the original; (b) is the filtered
image; (c) is the sorted coefficients; (d), (e) and (f) are the reconstructions
based on the largest 100, 300 and 500 coefficients, respectively. . . . . . . . 139
xxii
A.3 Edgelet transform of the wood grain image: (a) is the original; (d) is the
sorted coefficients; (b), (c), (e) and (f) are reconstructions based on the
largest 1 × 104 , 2 × 104 , 4 × 104 , 8 × 104 coefficients, respectively. . . . . . . 140
A.4 Edgelet transform of Lenna image: (a) is the original Lenna image; (b) is the
filtered version; (c) is the sorted largest 5, 000 coefficients out of 428032. (d),
(e) and (f) are the reconstructions based on the largest 1000, 2000 and 4000
coefficients, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
A.5 Vertices at scale j, for a 8 × 8 image with l = 1. The arrows shows the trend
of ordering. Integers outside are the labels of vertices. . . . . . . . . . . . . 144
A.6 Weights of pixels for one edgelet coefficient. . . . . . . . . . . . . . . . . . . 147
xxiii
xxiv
Nomenclature
Special sets
N ........................... Set of natural numbers {1, 2, 3, . . . }
R ........................... Set of real numbers
Z ........................... Set of integer numbers {. . . , −2, −1, 0, 1, 2, 3, . . . }
Operators
I ........................... Insert operator
S ........................... Downsample operator
Matrices
T ........................... Toeplitz matrix
H ........................... Hankel matrix
C ........................... Circulant matrix
Miscellaneous
f (·) . . . . . . . . . . . . . . . . . . . . . . . . . Function with domain R
f [·] . . . . . . . . . . . . . . . . . . . . . . . . . Function with domain Z
X ........................... Fourier transform of X
{· · ·} . . . . . . . . . . . . . . . . . . . . . . . . Ordered sequence
1 ........................... All one column vector
Pr. . . . . . . . . . . . . . . . . . . . . . . . . . . Probability
xxv
xxvi
List of Abbreviations
xxvii
LPF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Low-Pass Filter
LS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Least Square
LTI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Linear and Time Invariant
MAD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mean of the Absolute Deviation
MOF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Method of Frames
MOFDN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Method-of-Frames De-Noising
MP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Matching Pursuit
MPDN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Matching Pursuit De-Noising
MPEG . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Moving Picture Experts Group
MRA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Multiresolution Analysis
MRF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Markov Random Field
MSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mean-Square Error
OMP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Orthogonal Matching Pursuit
PCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Principal Component Analysis
RF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Random Field
RMS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Residual Mean Square
SAI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sparse Approximate Inverse
SIRP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Sparse Image Representation Pursuit
TFP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Time Frequency Plane
TFR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Time Frequency Representation
TVDN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Total Variation De-Noising
p.d.f. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . probability density function
xxviii
Chapter 1
Introduction
1.1 Overview
Recently, many new methods of image representation have been proposed, including wavelets,
cosine packets, brushlets, edgelets, ridgelets, and so on. Typically each of these is good for a
specific class of features, but not good for others. We propose a method of combining image
representations, more particularly, a method based on the 2-D wavelet transform and the
edgelet-like transform. The 2-D wavelet transform is good at capturing point singularities,
while the newly proposed edgelet-like transform is good at capturing linear singularities
(edges). Both transforms have fast algorithms for digital images.
To develop an efficient algorithm to solve this problem, we utilize insights from convex
optimization and the LSQR solver from computational linear algebra. These methods
combine to give a fast algorithm for our problem. Numerical experiments provide promising
results.
1
2 CHAPTER 1. INTRODUCTION
1.2 Outline
This thesis is organized as follows:
• Chapter 2 explains why sparse decomposition may lead to a more efficient image
coding and compression. It contains two parts: the first part gives a brief description
of the mathematical framework of a modern communication system; the second part
presents some quantitative results to explain (mainly in the asymptotic sense) why
sparsity in coefficients may lead to efficiency in image coding.
• Chapter 5 is a survey of iterative methods and explains why we choose LSQR. Some
alternative approaches are discussed.
• Chapter 6 presents some numerical simulations. They show that a combined approach
does provide a sparser decomposition than the existing approach that uses only one
transform.
• Appendix A documents the details about how to implement the exact edgelet trans-
form in a direct way. This algorithm has high complexity. Some examples are given.
This chapter explains the connection between sparse image decomposition and efficient
image coding. The technique we have developed promises to improve the efficiency of
transform coding, which is ubiquitous in the world of digital signal and image processing.
We start with an overview of image coding and emphasize the importance of transform
coding. Then we review Donoho’s work that answers in a deterministic way the question:
why does sparsity lead to good compression? Finally we review some important principles
in image coding.
Nowadays, it is almost impossible to review the field of theoretical signal processing without
mentioning the contributions of Shannon. In his 1948 paper [127], Shannon writes, “The
3
4 CHAPTER 2. SPARSITY IN IMAGE CODING
2. The transmitter, or encoder, transfers the message into a signal compatible with the
channel. In old-fashioned telephony, this operation consists of converting sound pres-
sure into a proportional electrical current. In telegraphy, we have an encoding system
to transfer a sequence of words to a sequence of dots, dashes and spaces. More complex
operations are applied to messages in modern communication systems.
3. The channel is the physical medium that conducts the transmitted signal. The math-
ematical abstraction of the transmission is a perturbation by noise. A typical as-
sumption on the channel is that the noise is additive, but this assumption can be
changed.
4. The receiver, or decoder, attempts to retrieve the message from the received signal.
Naturally, the decoding scheme depends on the encoding scheme.
Shannon’s model is both concrete and flexible. It has been proven efficient for model-
ing real world communication systems. From this model, some fundamental problems of
communication theory are apparent.
2.1. IMAGE CODING 5
INFORMATION
SOURCE TRANSMITTER CHANNEL RECEIVER DESTINATION
SIGNAL RECEIVED
SIGNAL
MESSAGE MESSAGE
NOISE
SOURCE
2. Given that the answer to the previous question is yes, how efficient can the system
(including both transmitter and receiver) be?
3. To achieve the maximum capacity, how should the transmitter and receiver be de-
signed?
To answer these three questions thoroughly would require a review of the entire field
of information theory. Here, we try to give a very brief, intuitive and non-mathematically-
rigorous description. A good starting point in information theory is the textbook [33] by
Cover and Thomas.
We need to introduce three concepts: entropy rate, channel capacity and distortion.
The entropy rate is a measure of complexity of the information source. Let us consider
only a discrete information source. At times n ∈ N, the information source emits a message
Xn from a discrete set X; thus, Xn ∈ X, for all n. So an infinite sequence of messages
X = {Xn , n ∈ N} is in X∞ . Suppose each Xn is drawn independently and identically
distributed from X with probability density function p(x) = Pr.{Xn = x}. Then the
entropy rate of this information source is defined as
H = H(X) = H(p(·)) = − p(x) log2 p(x). (2.1)
x∈X
6 CHAPTER 2. SPARSITY IN IMAGE CODING
Another intuitive way to view this is that if p(x) were a uniform distribution over X, then
the entropy rate H(p(·)) would be log2 |X|, where |X| is the cardinality (size) of the set X.
In (2.1) we can view 1/ log2 p(x) as an analogue to the cardinality such that the entropy is
the average logarithm of it.
The channel capacity, similarly, is the average of the logarithm of the number of bits
(binary numbers from {0, 1}) that can be reliably distinguished at receiver. Let Y denote
the output of the channel. The channel capacity is defined as
C = max I(X, Y ),
p(X)
where I(X, Y ) is the mutual information between channel input X and channel output Y :
The distortion is the statistical average of a (usually convex) function of the difference
between the estimation X̂ at the receiver and the original message X. Let d(x̂ − x) be the
measure of deviation. There are many options for function d: e.g., (1) d(x̂ − x) = (x̂ − x)2 ,
(2) d(x̂ − x) = |x̂ − x|. The distortion D is simply the statistical mean of d(x̂ − x):
D= d(x̂ − x)p(x)dx.
For (1), the distortion D is the residual mean square (RMS); for (2), the distortion D is
the mean of the absolute deviation (MAD).
The Fundamental Theorem for a discrete channel with noise [127] states that commu-
nication with an arbitrarily small probability of error is possible (from an asymptotic point
of view, by coding many messages at once) if H < C and impossible if H > C.
TRANSMITTER
SOURCE CHANNEL
CODER CODER
Shannon’s Fundamental Theorem of communication has two parts: the direct part and
the converse part. The direct part states that if the minimum achievable source coding
rate of a given source is strictly below the capacity of a channel, then the source can be
reliably transmitted through the channel by appropriate encoding and decoding operations;
the converse part states that if the source coding rate is strictly greater than the capacity of
channel, then a reliable transmission is impossible. Shannon’s theorem implies that reliable
transmission can be separated into two operations: source coding and channel coding. The
design of source encoding and decoding can be independent of the characteristic of the
channel; similarly, the design of channel encoding and decoding can be independent of the
characteristics of the source. Owing to the converse theorem, the reliable transmission is
doable either through a combination of separated source and channel coding or not possible
at all—whether it is a joint source channel coding or not. Note that in Shannon’s proof, the
Fundamental Theorem requires the condition that both source and channel are stationary
and memoryless. A detailed discussion about whether the Fundamental Theorem holds in
more general scenarios is given in [138].
One fact I must point out is that in statistical language, Shannon’s theorem is an
8 CHAPTER 2. SPARSITY IN IMAGE CODING
The purpose of source encoding is to transfer a message into a bit stream; and source
decoding is essentially an inverse operation. There is no need to overstate the ubiquity
of transform coding in modern digital communication systems. For example, JPEG is an
industry standard for still image compression and transmission. It implements a linear
transform—a two-dimensional discrete cosine transform (2-D DCT). The MPEG standard
is an international standard for video compression and transmission that incorporates the
JPEG standards for still image compression. The MPEG standard could be the standard
for future digital television.
Figure 2.3 gives a general depiction of a modern digital communication system with em-
bedded transform coding. The input message is x. The operator T is a transform operator.
Most of the time, T is a linear time invariant (LTI) transform. After the transform, we
get coefficients y. Operator Q is a quantizer—it transfer continuous variables (in this case,
they are coefficients after transform T ) into discrete values. Without loss of generality, we
can assume each message is transferred into a bit stream of finite length. The operator E
is an entropy coder. It further shortens the bit stream by using, for example, run-length
coding for pictures [23], Huffman coding, arithmetic coding, or Lempel-Ziv coding, etc. Op-
erators T , Q and E together make the encoder. We assume that the bit stream is perfectly
transmitted through an error-corrected channel. The operators E −1 and Q−1 are inverse
operators of the entropy coder E and the quantizer Q, respectively. The operator U is the
reconstruction transform of transform T . If transform T is invertible, then the transform
U should be the inverse transform of T , or equivalently, U = T −1 . The final output at the
receiver is the estimation, denoted by x̂, of the original message. Operators E −1 , Q−1 and
U together form the decoder.
2.1. IMAGE CODING 9
ERROR-CORRECTED
CHANNEL
x- y- - E - E −1 - Q−1 - U
x̂ -
T Q
Figure 2.3: Structure of a communication system with a transform coder. The signal is
x. Symbol T denotes a transform operator. The output coefficients vector (sequence) is
y. Symbols Q and E denote a quantizer and an entropy coder, respectively. Symbols Q−1
and E −1 denote the inverse operators (if they exist) of operations Q and E. Symbol U
stands for the reconstruction transform of transform T . Symbol x̂ stands for the estimated
(/received/recovered) signal.
The idea of transform coding is simply to add a transform at the sender before a quantizer
and to add the corresponding reconstruction (or sometimes inverse) transform at the receiver
after an inverse of the quantizer. Transform coding improves the overall efficiency of a
communication system if the coefficient vector in the transform domain (space made by
coefficient vectors) is more sparse than the signal in the signal domain (space made by
original signals). In general, the sparser a vector is, the easier it is to compress and quantize.
So when the coefficient vector is sparser, it is relatively easy to quantize and compress in the
transform domain than in the signal domain. The transform should have a corresponding
reconstruction transform, or even be invertible, as we should be able to reconstruct the signal
at the receiver. The efficiency of the communication system in Figure 2.3 is determined
mostly by the effectiveness of the transform coder.
Now we talk about the history of transform coding in digital communication. The orig-
inal idea was to use transforms to decorrelate dependent signals; later some researchers
successfully explored the possibility of using transforms to improve the efficiency of com-
munication systems. The next paragraph describes this history.
In 1933, in the Journal of Educational Psychology, Hotelling [82] first presented a decor-
relation method for discrete data. This is the starting point of a popular methodology
now called principal component analysis (PCA) in the statistics community. The analogous
transform for continuous data was obtained by Karhunen in 1947 and Loéve in 1948. The
two papers by Karhunen and Loéve are not in English, so they are rarely cited. Readers can
find them in the reference in [126]. In 1956, Kramer and Mathews [115] introduced a method
for transmitting correlated, continuous-time, continuous-amplitude signals. They showed
that the total bandwidth (which in theory is equivalent to the capacity of the channel)
10 CHAPTER 2. SPARSITY IN IMAGE CODING
necessary to transmit a set of signals with prescribed fidelity can be reduced by transmit-
ting a set of linear combinations of the signals, instead of the original signals. Assuming
Gaussian signals and RMS or mean-squared error (MSE) distortion, the Karhunen-Loéve
transform (KLT) is optimal. Later, Huang and Schultheiss [84, 85] extended this idea to
a system that includes a quantizer. Nowadays, PCA has become a ubiquitous technique
in multivariate analysis. Its equivalent—KLT—has become a well-known method in the
electrical engineering community.
To find a transform that will generate even sparser coefficients than most of the existing
transform is the ultimate goal of this thesis. Our key idea is to combine several state-of-the-
art transforms. The combination will introduce redundancy, but at the same time it will
provide flexibility and potentially lead to a decomposition that is sparser than one with a
single transform.
When describing image coding, we cannot avoid mentioning three topics: quantization,
entropy coding and predictive coding. Each of the following subsubsections is dedicated to
one of these topics.
Quantization
All quantization schemes can be divided into two categories: scalar quantization and vector
quantization. We can think of vector quantization as an extension of scalar quantization in
high-dimensional space. Gray and Neuhoff 1998 [76] provides a good review paper about
quantization. Two well-written textbooks are Gersho and Gray 1992 [70] and Sayood 1996
[126].
The history of quantization started from Pulse-code modulation (PCM), which was
patented in 1938 by Reeves [122]. PCM was the first digital technique for conveying an
analog signal over an analog channel. Twenty-five years later, a noteworthy article [37] was
written by Reeves and Deloraine. It gave a historical review and an appraisal of the future
of PCM. Some predictions in this article turned out to be prophetic, especially the one
about the booming of digital technology.
Here we give a brief description of scalar quantization. A scalar quantization is a map-
ping from a real-valued source to a discrete subset; symbolically,
Q : R → C,
2.1. IMAGE CODING 11
where Q stands for the quantization operator and set C is a discrete subset of the real
number set R:
The procedure of designing a quantizer can be divided into two steps: partitioning
and representers selection. In the partitioning step, the real number set R is partitioned
into a group of subsets c1 , c2 , . . . , cK —Note that this K is the same K as in (2.2)—so
that these subsets are mutually exclusive (ci ∩ cj = ∅, ∀i = j) and also comprehensive
(∪i=1,2,... ,K ci = R). Each of these ci ’s are called a cell. In the step of choosing representers,
each cell ci is associated with a real value yi (the representer), such that when x ∈ ci , the
quantization of x (denoted by q(x)) is equal to yi . Equivalently,
x ∈ ci =⇒ q(x) = yi , i = 1, 2, . . . , K.
In the measurement of quantization efficiency, two quantities, distortion and rate, are
important. The difference between a quantizer input x and a quantizer output q(x), namely,
q(x) − x, is called the quantization error. The distortion of a quantizer is defined as a
statistical average of a non-negative convex function, say d, of the quantization error:
The rate, in the finite number of cells case, is the logarithm of the number of cells. For
example, in (2.2) the rate is log(K). In a given system, the rate and distortion are two
competing quantities. The larger the rate is, the smaller the minimum distortion could be.
And vice versa.
Consider the MSE distortion measure, D(q) = E[(q(x) − x)2 ]. For each step in quanti-
zation, there is a necessary condition for the optimality of the design. So we have in total
two conditions:
• Nearest neighbor classification: This is a rule for partitioning. Consider the codebook
to be fixed. In encoding a source sample x, one should choose the codebook element
12 CHAPTER 2. SPARSITY IN IMAGE CODING
q(x) = argmin |x − yi |.
yi ∈ C
yi = E[x|x ∈ ci ], i = 1, 2, . . . , K.
Given a source with a known statistical probability density function (p.d.f.), designing a
quantizer that satisfies both of the optimality conditions is not a trivial task. In general, it
is done via an iterative scheme. A quantizer designed in this manner is called a Lloyd-Max
quantizer [94, 104].
A uniform quantizer simply takes each cell ci as an interval of equal length; for example,
for fixed ∆, ci is (i∆ − ∆
2 , i∆ + ∆
2 ]. (Note that the number of cells here is not finite.) The
uniform quantizer is easy to implement and design, but in general not optimal.
A detailed discussion about quantization theory is easy to find and beyond the scope of
this thesis. We skip it.
Entropy Coding
Entropy coding means compressing a sequence of discrete values from a finite set into a
shorter sequence combined from elements of the same set. Entropy coding is a lossless
coding, which means that from the compressed sequence one can perfectly reconstruct the
original sequence.
There are two important approaches in entropy coding. One is a statistical approach.
This category includes Huffman coding and arithmetic coding. The other is a dictionary
based approach. The Lempel-Ziv coding belongs to this category. For detailed descriptions
of these coding schemes we refer to some textbooks, notably Cover and Thomas 1991 [33],
Gersho and Gray 1992 [70] and Sayood 1996 [126].
2.2. SPARSITY AND COMPRESSION 13
Predictive Coding
For an image transmission system applying transform coding, the general scheme is depicted
in Figure 2.4. First of all, a transform is applied to the image, and a coefficient vector y is
obtained. Generally speaking, the sparser the vector y is, the easier it can be compressed.
The coefficient vector y is quantized, compressed, and transmitted to the other end of
the communication system. At the receiver, the observed vector ỹ should be a very good
approximation to the original coefficient vector y. An inverse transform is applied to recover
the image.
We have not answered the question on how to quantify the sparsity of a vector and why
the sparsity leads to a “good” compression. In the next section, to answer these questions,
we introduce the work of Donoho that gives a quantification of sparsity and its implication
in effective bit-level compression.
IMAGE
?
TRANSFORM
?
QUANTIZATION/COMPRESSION
BITS
?
CHANNEL
BITS
?
DECOMPRESSION/RECONSTRUCTION
ỹ
?
INVERSE TRANSFORM
IMAGE
?
Figure 2.4: Diagram of image transmission using transform coding. The sequence (or
vector) y is the coefficient sequence (or vector) and ỹ is the recovered coefficient sequence
(or vector).
2.2. SPARSITY AND COMPRESSION 15
1. What is the weak lp norm and what are the scientific reasons to invent such a norm?
(Answered in Section 2.2.1.)
We start with the definition of the weak lp norm, for 1 ≤ p ≤ 2, then describe its connec-
tion with compression number, measure of numerosity and rate of recovery, and finally we
introduce the concept of critical index, which plays an important role in the next section
(Section 2.2.2).
16 CHAPTER 2. SPARSITY IN IMAGE CODING
Note that this is a quasi-norm. (A norm satisfies a triangular inequality, x+y ≤ x+y;
a quasi-norm satisfies a quasi-triangular inequality, x+y ≤ K(x+y) for some K > 1.)
The weak lp norm has a close connection with three other measures that are related to vector
sparsity.
Numerosity
A straightforward way to measure the sparsity of a vector is via its numerosity: the number
of elements whose amplitudes are above a given threshold δ. In a more mathematical
language, for a fixed real value δ, the numerosity is equal to #{i : |θi | > δ}. The following
lemma is cited from [47].
From the above lemma, a small weak lp norm leads to a small number of elements that
are significantly above zero. Since numerosity, which basically counts significantly large
elements, is an obvious way to measure the sparsity of a vector, the weak lp norm is a
measure of sparsity.
Compression Number
Again, the compression number is based on the sorted amplitudes. In an orthogonal basis,
we have isometry. If we perform a thresholding scheme by keeping the coefficients associated
2.2. SPARSITY AND COMPRESSION 17
with the largest n amplitudes, then the compression number c(n) is the square root of the
RSS distortion of the signal reconstructed by the coefficients with the n largest amplitudes.
The following result can be found in [47].
Lemma 2.2 For any sequence θ, if m = 1/p − 1/2, the following inequality is true:
c(N ) ≤ αp N −m |θ|wlp , N ≥ 1,
From the above lemma, a small weak lp norm implies a small compression number.
Rate of Recovery
The rate of recovery comes from statistics, particularly in density estimation. For a sequence
θ, the rate of recovery is defined as
∞
r() = min{θi2 , 2 }.
i=1
Lemma 2.3 For any sequence θ, if r = 1 − p/2, the following inequality is true:
This implies that a small weak lp norm leads to a small rate of recovery. In some cases
(for example, in density estimation) we choose rate of recovery as a measure of sparsity.
The weak lp norm is therefore a good measure of sparsity too.
Lemma 1 in [46] shows that all these measures are equivalent in an asymptotic sense.
Critical Index
In order to define the critical index of a functional space, we need to introduce some new
notation. A detailed discussion of this can be found in [47]. Suppose Θ is the functional
space that we are considering. (In the transform coding scenario, the functional space Θ
includes all the coefficient vectors.) An infinite-length sequence θ = {θi : i ∈ N} is in a weak
18 CHAPTER 2. SPARSITY IN IMAGE CODING
lp space if and only if θ has a finite weak lp norm. For fixed p, the statement “Θ ⊂ wlp ”
simply means that every sequence in Θ has a finite weak lp norm. The critical index of a
functional space Θ, denoted by p∗ (Θ), is defined as the infimum of p such that the weak lp
space includes Θ:
In Section 2.2.2, we describe the linkage between critical index and optimal exponent,
where the latter is a measure of efficiency of “optimal” coding in a certain functional space.
Before we move into the asymptotic discussion, the following result provides insight into
coding a finite-length sequence. The key idea is that the faster the finite sequence decays,
the fewer bits required to code the sequence. Here the sparsity is measured by the decay of
the sorted amplitudes, which in the 1-D case is the decay of the sorted absolute values of a
finite-length sequence.
In general, it is hard to define a “universal” measure of speed of decay for a finite-length
sequence. We consider a special case. The idea is depicted in Figure 2.5. Suppose that
two curves of sorted amplitudes of two sequences (indicated by A and B in the figure) have
only one intersection. We consider the curve that is lower at the tail part (indicated by
B in the figure) corresponds to a sparser sequence. (Obviously, this is a simplified version
because there could be more than one intersection.) We prove that by using the most
straightforward coding scheme, which is to take an identical uniform scalar quantizer at
every coordinate, we need fewer bits to code B than to code A.
To be more specific, without loss of generality, suppose we have two normalized non-
increasing sequences. Both of these sequences have L2 norm equal to 1. One decays “faster”
than the other, in the sense that when after a certain index, the former sorted-amplitude
sequence is always below the latter one, as in Figure 2.5. (Curve A is always above curve
B after the intersection.) Suppose we deploy a uniform scalar quantizer with the same
quantization parameter q on every coordinate. Then the number of bits required to code
the ith element θi in the sequence θ is no more than log2 [θi /q] + 1 (where [x] is the closest
integral value to x). Hence the total number of bits to code vector θ = {θi : 1 ≤ i ≤ N }
is upper bounded by N i=1 log2 [θi /q] + N . Since q and N are constants, we only need to
2.2. SPARSITY AND COMPRESSION 19
log(|c| )
(i)
A
Figure 2.5: Intuition to compare the speeds of decay for two sequences. A sequence whose
sorted-amplitudes-curve B is, as we consider, sparser than the sequence whose sorted-
amplitudes-curve is A.
N
consider i=1 log |θi |. The following theorem shows that a “sparser” sequence needs fewer
bits to code.
First suppose {xi : 1 ≤ i ≤ N } and {yi : 1 ≤ i ≤ N } are two non-increasing positive
real-valued sequences with a fixed l2 norm:
N
x1 ≥ x2 ≥ . . . ≥ xN > 0, 2
i=1 xi = C;
N
y1 ≥ y2 ≥ . . . ≥ yN > 0, 2
i=1 yi = C;
Theorem 2.1 If sequence x majorizes sequence y and both of them have the same 2 norm,
N
then Ni=1 log xi ≥ i=1 log yi .
In the previous section, we described the weak lp norm and critical index. It is interesting
to note that the critical index is closely related to the asymptotics of compression in a
functional space. To describe this relationship, we need to introduce a new concept: optimal
20 CHAPTER 2. SPARSITY IN IMAGE CODING
exponent. We first introduce the definition of optimal exponent, then cite the result that
states the connection between critical index and optimal exponent.
Optimal Exponent
The -entropy of Θ, H (Θ), is to within one bit the minimum number of bits required
to code in space Θ with distortion less than . If we apply the nearest-neighbor coding,
the total number of possible cases needed to record is N (, Θ), which should take no more
than log2 N (, Θ) bits. (The value x is the smallest integer that is no smaller than x.)
Obviously the -entropy of Θ, which is denoted by H (Θ), is within 1 bit of the total number
of bits required to code in functional space Θ.
It is not hard to observe that when → 0, H (Θ) → +∞. Now the question is how fast
the -entropy H (Θ) increases when goes to zero. The optimal exponent defined in [47] is
a measure of the speed of increment:
We refer readers to the original paper [47] for a more detailed discussion.
Asymptotic Result
The key idea is that there is a direct link between the critical index p∗ and the optimal
exponent α∗ . The following is cited from [47], Theorem 2.
Theorem 2.2 Let Θ be a bounded subset of l2 that is solid, orthosymmetric, and minimally
2.3. DISCUSSION 21
Moreover, coder-decoder pairs achieving the optimal exponent of code length can be derived
from simple uniform quantization of the coefficients (θi ), followed by simple run-length
coding.
2.2.3 Summary
1. Critical index measures the sparsity of a sequence space, which usually is a subspace
of l2 . In an asymptotic sense, the smaller the critical index is, the faster the sequence
decays; and, hence, the sparser the sequence is.
2. Optimal exponent measures the efficiency of the best possible coder in a sequence
space. The larger the optimal exponent is, the fewer the bits required for coding in
this sequence space.
3. The optimal exponent and the critical index have an equality relationship that is
described in (2.4). When the critical index is small, the optimal exponent is large.
Together with the previous two results, we can draw the main conclusion: the sparser
the sequences are in the sequence space, the fewer the bits required to code in the
same space.
2.3 Discussion
Based on the previous discussion, in an asymptotic sense, instead of considering how many
bits are required in coding, it is equivalent to study the sparsity of the coefficient sequences.
From now on, for simplicity, we only consider the sparsity of coefficients. Note that the
sparsity of coefficients can be empirically measured by the decay of sorted amplitudes. For
infinite sequences, the sparsity can be measured by the weak lp norm.
We required the decoding scheme to be computationally cheap, while the encoding
scheme can be computationally expensive. The reason is that in practice (for example, in
22 CHAPTER 2. SPARSITY IN IMAGE CODING
the broadcasting business), the sender (TV station) can usually afford expensive equipment
and a long processing time, while the receiver (single television set) must have cheap and
fast (even real-time) algorithms.
Our objective is to find sparse decompositions. Following the above rule, in our algo-
rithm we will tolerate high complexity in decomposing, but superposition must be fast and
cheap. Ideally, the order of complexity of superpositioning must be no higher than the order
of complexity of a Fast Fourier Transform (FFT). We choose FFT for comparison because
it is a well-known technique and a milestone in the development of signal processing. Also
by ignoring the logarithm factor, the order of complexity of FFT is almost equal to the
order of complexity of copying a signal or image from one disk to another. In general, we
can hardly imagine any processing scheme that can have a lower order of complexity than
just copying. We will see that the order of complexity of our superposition algorithm is
indeed no higher than the order of complexity of doing an FFT.
2.4 Proof
N
N
log x2i ≥ log yi2 . (2.5)
i=1 i=1
N
f (x) = log x2i .
i=1
(n−1)
l = min{i : yi > xi }.
2.4. PROOF 23
The index l does not exist if and only if the two sequences {xi , 1 ≤ i ≤ N } and
{yi , 1 ≤ i ≤ N } are the same.
(b) Pick an index u such that
(n−1)
u = max{i : yi < xi }.
Again, the index u does not exist if and only if the two sequences {xi , 1 ≤ i ≤ N }
and {yi , 1 ≤ i ≤ N } are the same.
(c) Take
(n−1) 2
ε = min{yl2 − (xl ) , (x(n−1)
u )2 − yu2 }.
2. The difference u − l has an integral value bounded between 1 and N − 1. At each step,
the difference u − l should reduce by at least one. Based on this, the above procedure
can repeat at most N − 1 times.
3. It is clear that f (x(n−1) ) > f (x(n) ), n ≥ 1. And because f (x) = f (x(0) ) and f (y) is
equal to the last function value in f (x(n) ), n ≥ 1, we have f (x) ≥ f (y), which is the
inequality (2.5).
This chapter is a summary of early achievements and recent advances in the design of
transforms. It also describes features that these transforms are good at processing.
Some mathematical slogans are given, together with their quantitative explanations. As
many of these topics are too large to be covered in a single chapter, we give pointers in the
literature. The purpose of this chapter is to build a concrete foundation for the remainder
of this thesis.
One particular message that we want to deliver is: “Each transform is good for one
particular phenomenon, but not good for some others; at the same time, a typical image is
made by a variety of phenomena. It is natural to think about how to combine different trans-
forms, so that the combined method will have advantages from each of these transforms.”
This belief is the main motivation for this thesis.
We give emphasis to discrete algorithms and, moreover, fast linear algebra, because the
possibility of efficient implementation is one of our most important objectives.
Note that there are two important ideas in signal transform methods: decomposition and
distribution. The idea of decomposition, also called atomic decomposition, is to write the
image/signal as a superposition (which is equivalently a linear combination) of pre-selected
atoms. These atoms make a dictionary, and the dictionary has some special property;
for example, orthonormality when a dictionary is made by a single orthonormal basis, or
tightness when a dictionary is made by a single tight frame. The second idea is to map the
signal into a distribution space. An example is the Wigner-Ville transform, which maps a
25
26 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
1-D signal into a distribution on the time-frequency plane (which is 2-D). Since this thesis
mainly emphasizes decomposition, from now on “transform” means only the one in the
sense of decomposition.
Some statements could be seemingly subjective. Given the same fact (in our case,
for example, it could be a picture), different viewers may draw different conclusions. We
manage to be consistent with the majority’s point of view. (Or at least we try to.)
Key message
Before we move into the detailed discussion, we would like to summarizes the key results.
Readers will see that they play an important role in the following chapters. Table 3.1
summarize the correspondence between three transforms, three image features, and orders
of complexity of their fast algorithms. Note that the 2-D DCT is good for an image with
homogeneous components. For an N by N image, its fast algorithm has O(N 2 log N )
order of complexity. The two-dimensional wavelet transform is good for images with point
singularities. For an N by N image, its fast algorithm has O(N 2 ) order of complexity. The
edgelet transform is good for an image with line singularities. For an N by N image, its
order of complexity is O(N 2 log N ), which is the same as the order of complexity for the
2-D DCT.
Table 3.1: Comparison of transforms and the image features that they are good at process-
ing. The third column lists the order of computational complexity for their discrete fast
algorithms.
In the main body of this chapter, we show many figures. The most important reason is that
this project is an image processing project, and showing figures is the most intuitive way
to illustrate points.
27
Some of the figures show the basis functions for various transforms. The reason for
illustrating these basis functions is the following: the ith coefficient of a linear transform is
the inner product of the ith basis function and the image:
where ·, · is the inner product of two functions. Hence the pattern of basis functions
determines the characteristics of the linear transform.
In functional analysis, a basis function is called a Riesz representer of the linear trans-
form.
Notations
We follow some conventions in signal processing and statistics. A function X(t) is a contin-
uous function with a continuous time variable t. A function X[k] is a continuous function
is the Fourier transform of X.
with variable k that only takes integral values. Function X
We use ω to denote a continuous frequency variable.
More specifics of notation can be found in the main body of this chapter.
Organization
The rest of this chapter is organized as follows. The first three sections describe the three
most dominant transforms that are used in this thesis: Section 3.1 is about the discrete
cosine transform (DCT); Section 3.2 is about the wavelet transform; and Section 3.3 is
about the edgelet transform. Section 3.4 is an attempt at a survey of other activities.1
Section 3.5 contains some discussion. Section 3.6 gives conclusions. Section 3.7 contains
some proofs.
1
It is interesting to note that there are many other activities in this field. Some examples are the Gabor
transform, wavelet packets, cosine packets, brushlets, ridgelets, wedgelets, chirplets, and so on. Unfortu-
nately, owing to space constraints, we can hardly provide many details. We try to maintain most of the
description at an introductory level, unless we believe our detailed description either gives a new and useful
perspective or is essential for some later discussions.
28 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
These three items comprise the content of Section 3.1.1. Section 3.1.2 switches to the
discrete cosine transform (DCT). It has two ingredients:
1. definition,
2. fast algorithms.
For fast algorithms, the theory follows two lines: (1) An easy way to implement fast DCT is
to utilize fast DFT; the fast DFT is also called the fast Fourier transform (FFT). (2) We can
usually do better by considering the sparse matrix factorization of the DCT matrix itself.
(This gives a second method for finding fast DCT algorithms.) The algorithms from the
second idea are computationally faster than those from the first idea. We summarize recent
advances for fast DCT in Section 3.1.2. Section 3.1.3 gives the definitions of Discrete Sine
Transform (DST). Section 3.1.4 provides a framework for homogeneous images, and explains
when a transform is good at processing homogeneous images. The key idea is that the
DCT almost diagonalizes the covariance matrix of certain Gaussian Markov random fields.
Section 3.1.4 shows that the 2-D DCT is a good transform for images with homogeneous
components.
Definition
In order to introduce the DFT, we need to first introduce the continuous Fourier transform.
From the book by Y. Meyer [107, Page 14], we learn that the idea of continuous Fourier
transform dates to 1807, the year Joseph Fourier asserted that any 2π-periodic function
in L2 , which is a functional space made by functions whose square integral is finite, is
3.1. DCT AND HOMOGENEOUS COMPONENTS 29
a superposition of sinusoid functions. Nowadays, Fourier analysis has become the most
powerful tool in signal processing. Every researcher in science and engineering should know
it.
Since digital signal processing (DSP) has become the mainstream in signal processing,
and only discrete transforms can be implemented in the digital domain, it is more interesting
to consider the discrete version of Fourier transform.
For a finite sequence X[0], X[1], . . . , X[N − 1], where N is a given integer, the discrete
Fourier transform is defined as
N −1
= √1
X[l] X[k]e−i N kl ,
2π
0 ≤ l ≤ N − 1; (3.1)
N k=0
and consequently
N −1
1 i 2π kl
X[k] = √ X[l]e N , 0 ≤ k ≤ N − 1.
N l=0
l=0,1,2,... ,N −1
= DFTN · X.
X
are N -
Apparently DFTN is an N by N square symmetric matrix. Both X and X
dimensional vectors.
It is well known that the matrix DFTN is orthogonal, which is equivalent to saying
that the inverse of DFTN is its complex conjugate transpose. Thus, DFT is an orthogonal
transform and consequently it is an isometric transform. This result is generally attributed
to Parseval.
Since the concise statements and proofs can be found in many textbooks, we skip the
proof.
The preponderant reason that Fourier analysis is so powerful in analyzing linear time invari-
ant system and cyclic-stationary time series is the fact that Fourier series are the eigenvectors
of matrices associated with these transforms. These transforms are ubiquitous in DSP and
time series analysis.
Before we state the key result, let us first introduce some mathematical notation:
Toeplitz matrix, Hankel matrix and circulant matrix.
Suppose AN ×N is an N by N real-valued matrix: A = {aij }1≤i≤N,1≤j≤N with all aij ∈ R.
The matrix A is a Toeplitz matrix when the elements on the diagonal, and the rows that
are parallel to the diagonal, are the same, or mathematically aij = ti−j [74, page 193]. The
following is a Toeplitz matrix:
⎛ ⎞
t0 t−1 ... t−(n−2) t−(n−1)
⎜ ⎟
⎜ t ⎟
⎜ 1 t0 . . . t−(n−3) t−(n−2) ⎟
⎜ . .. .. .. ⎟
T =⎜
⎜ .
. .
..
. . . ⎟
⎟ .
⎜ ⎟
⎜ tn−2 tn−3 . . . t0 t−1 ⎟
⎝ ⎠
tn−1 tn−2 . . . t1 t0
N ×N
3.1. DCT AND HOMOGENEOUS COMPONENTS 31
A matrix A is a Hankel matrix when its elements satisfying aij = hi+j−1 . A Hankel
matrix is symbolically a flipped (left-right) version of a Toeplitz matrix. The following is a
Hankel matrix:
⎛ ⎞
h1 h2 ... hn−1 hn
⎜ ⎟
⎜ h ⎟
⎜ 2 h3 . . . hn hn+1 ⎟
⎜ . .. .. .. ⎟
H=⎜
⎜ .
. .
..
. . . ⎟
⎟ .
⎜ ⎟
⎜ hn−1 hn . . . h2n−3 h2n−2 ⎟
⎝ ⎠
hn hn+1 . . . h2n−2 h2n−1
N ×N
A matrix A is called a circulant matrix [74, page 201] when aij = c{i−j mod N } , where
x mod N is the non-negative remainder after a modular division. The following is a circulant
matrix:
⎛ ⎞
c0 cN −1 . . . c2 c1
⎜ ⎟
⎜ c ⎟
⎜ 1 c0 ... c3 c2 ⎟
⎜ . .. .. .. ⎟
C=⎜
⎜ .
. . ... . . ⎟
⎟ . (3.2)
⎜ ⎟
⎜ cN −2 cN −3 ... c0 cN −1 ⎟
⎝ ⎠
cN −1 cN −2 ... c1 c0
N ×N
Theorem 3.1 For an N × N circulant matrix C as defined in (3.2), the Fourier series
{e−i N kl : k = 0, 1, 2, . . . , N − 1}, for l = 0, 1, 2, . . . , N − 1, are eigenvectors of the ma-
2π
√ √ √
trix C, and the Fourier transforms of the sequence { N c0 , N c1 , . . . , N cN −1 } are its
eigenvalues.
1. For a cyclic linear time-invariant (LTI) system, the transform matrix is a circulant
matrix. Thus, Fourier series are the eigenvectors of a cyclic linear time invariant
system, and the Fourier transform of the square root N amplified impulse response is
32 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
the set of eigenvalues of this system. Note: when the signal is infinitely long and the
impulse response is relatively short, we can drop the cyclic constraint. Since an LTI
system is a widely assumed model in signal processing, the Fourier transform plays
an important role in analyzing this type of system.
From the above two examples, we see that FT plays an important role in both DSP and
time series analysis.
The fast Fourier transform (FFT) is one of the truly great computational devel-
opments of this century. It has changed the face of science and engineering so
2
No long-range dependency means that if two locations are far away in the series, then the corresponding
two random variables are nearly independent.
3.1. DCT AND HOMOGENEOUS COMPONENTS 33
There are many ways to derive the FFT. The following theorem, which is sometimes
called Butterfly Theorem, seems the most understandable way to describe why FFT works.
The author would like to thank Charles Chui [30] for being the first fellow to introduce this
theorem to me. The original idea started from Cooley and Tukey [32].
To state the theorem, let us first give some notation for matrices. Let P1 and P2 denote
two permutation matrices satisfying for any mn-dimensional vector x = (x0 , x1 , . . . , xmn−1 )
∈ Rmn , we have,
⎛ ⎛ ⎞ ⎞
x0
⎜ ⎜ ⎟ ⎟
⎜ ⎜ xn ⎟ ⎟
⎜ ⎜ ⎟ ⎟
⎜ ⎜ .. ⎟ ⎟
⎜ ⎜ ⎟ ⎟
⎜ ⎝ . ⎠ ⎟
⎜ x(m−1)n ⎟
⎜ ⎛ ⎞
m×1 ⎟
⎜ ⎟
⎜ x1 ⎟
⎜ ⎜ ⎟ ⎟
⎛ ⎞ ⎜ ⎜ x1+n ⎟ ⎟
⎜ ⎜
⎜
⎟
⎟
⎟
x0 ⎜ ⎜ .. ⎟ ⎟
⎜ ⎟ ⎜ ⎝ . ⎠ ⎟
⎜ x1 ⎟ ⎜ ⎟
P1 ⎜
⎜ .. ⎟
⎟ =⎜
⎜ x1+(m−1)n ⎟;
⎟ (3.3)
⎝ . ⎠ ⎜ m×1 ⎟
⎜ • ⎟
xmn−1 ⎜ ⎟
mn×1 ⎜ • ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎛
•
⎞
⎟
⎜ ⎟
⎜ x(n−1) ⎟
⎜ ⎜ ⎟ ⎟
⎜ ⎜ ⎟ ⎟
⎜ ⎜ x(n−1)+n ⎟ ⎟
⎜ ⎜ ⎟ ⎟
⎜ ⎜ .. ⎟ ⎟
⎝ ⎝ . ⎠ ⎠
x(n−1)+(m−1)n
m×1
34 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
and
⎛ ⎛ ⎞ ⎞
x0
⎜ ⎜ ⎟ ⎟
⎜ ⎜ xm ⎟ ⎟
⎜ ⎜ ⎟ ⎟
⎜ ⎜ .. ⎟ ⎟
⎜ ⎜ ⎟ ⎟
⎜ ⎝ . ⎠ ⎟
⎜ x(n−1)m ⎟
⎜ ⎛ ⎞
n×1 ⎟
⎜ ⎟
⎜ x1 ⎟
⎜ ⎜ ⎟ ⎟
⎛ ⎞ ⎜ ⎜ x1+m ⎟ ⎟
⎜ ⎜
⎜
⎟
⎟
⎟
x0 ⎜ ⎜ .. ⎟ ⎟
⎜ ⎟ ⎜ ⎝ . ⎠ ⎟
⎜ x1 ⎟ ⎜ ⎟
P2 ⎜
⎜ .. ⎟
⎟ =⎜
⎜ x1+(n−1)m ⎟.
⎟ (3.4)
⎝ . ⎠ ⎜ n×1 ⎟
⎜ • ⎟
xmn−1 ⎜ ⎟
mn×1 ⎜ • ⎟
⎜ ⎟
⎜ ⎟
⎜ ⎛
•
⎞
⎟
⎜ ⎟
⎜ x(m−1) ⎟
⎜ ⎜ ⎟ ⎟
⎜ ⎜ ⎟ ⎟
⎜ ⎜ x(m−1)+m ⎟ ⎟
⎜ ⎜ ⎟ ⎟
⎜ ⎜ .. ⎟ ⎟
⎝ ⎝ . ⎠ ⎠
x(m−1)+(n−1)m
n×1
Obviously, we have P1 P1T = Imn and P2 P2T = Imn , where Imn is an mn × mn identity
matrix. Like all permutation matrices, P1 and P2 are orthogonal matrices with only one
nonzero element in each row and column, and these nonzero elements are equal to one.
Let Ωmn denote a diagonal matrix, and let In be an n × n identity matrix. For ωmn =
−i mn
2π
e , the matrix Ωmn is
⎡ ⎤
In
⎢ ⎛ ⎞ ⎥
⎢ 1 ⎥
⎢ ⎜
⎜
⎟
⎟ ⎥
⎢ ⎜
1×1
ωmn ⎟ ⎥
⎢ ⎜ ⎟ ⎥
⎢ ⎜ .. ⎟ ⎥
⎢ ⎝ . ⎠ ⎥
⎢ ⎥
⎢ ωmn
(n−1)×1
⎥
Ωmn =⎢
⎢ ..
⎥ . (3.5)
⎥
⎢ . ⎥
⎢ ⎛ ⎞ ⎥
⎢ 1 ⎥
⎢ ⎜ ⎟ ⎥
⎢ ⎜ ωmn
1×(m−1)
⎟ ⎥
⎢ ⎜
⎜
⎟
⎟
⎥
⎢ ⎜ .. ⎟ ⎥
⎣ ⎝ . ⎠ ⎦
(n−1)×(m−1)
ωmn
3.1. DCT AND HOMOGENEOUS COMPONENTS 35
Armed with these definitions, we are now able to state the Butterfly Theorem. The
following is a formal statement.
Theorem 3.2 (Butterfly Theorem) Consider any two positive integers m and n. Let
DFTmn , DFTm and DFTn be the mn-, m- and n- point discrete Fourier transform ma-
trices respectively. The matrices P1 , P2 and Ωmn are defined in (3.3), (3.4) and (3.5),
respectively. The following equality is true:
⎡ ⎤ ⎡ ⎤
⎢ ⎥ ⎢ ⎥
⎢ DFTm ⎥ ⎢ DFTn ⎥
⎢ ⎥ ⎢ ⎥
⎢ .. ⎥ ⎢ .. ⎥ T
DFTmn = P1T ⎢ . ⎥ P2 Ωmn ⎢ . ⎥ P2 . (3.6)
⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢ ⎥
⎣ DFTm ⎦ ⎣ DFTn ⎦
n m
The proof would be lengthy and we omit it. At the same time, it is not difficult to find
a standard proof in the literature; see, for example, [136] and [134].
Why does Theorem 3.2 lead to a fast DFT? To explain this, we first use (3.6) in the
following way. Suppose N is even. Letting m = N/2 and n = 2, we have
⎡ ⎤
⎡ ⎤
⎢ ⎥
⎢ DFT ⎥
⎢ ⎥ ⎢ 2 ⎥
⎢
T ⎢ DFT ⎥ ⎢ ⎥ T
DFTN = P1 ⎢
N/2 ⎥ P2 ΩN ⎢ .. ⎥ P2 . (3.7)
⎥ ⎢ . ⎥
⎣ DFTN/2 ⎦ ⎢ ⎥
⎢ DFT2 ⎥
⎣ ⎦
2
N/2
Let C(N ) denote the number of operations needed to multiply a vector with matrix DFTN .
Operations include shifting, scalar multiplication and scalar addition. Since P1 and P2 are
permutation matrices, it takes N shifting operations to do vector-matrix-multiplication with
P1T , P2 and P2T . Note that DFT2 is a 2 × 2 matrix. Obviously, it takes 3N operations (2N
36 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
Finally, since ΩN is a diagonal matrix, it takes N multiplications to multiply with it. From
all the above, together with (3.7), we can derive the following recursive relationship:
C(N ) = N + 2C(N/2) + N + N + 3N + N
= 2C(N/2) + 7N. (3.8)
If N = 2m , we have
from (3.8)
C(N ) = 2C(N/2) + 7N
= ...
= 2m−1 C(2) + 7(m − 1)N
= 7N log2 N − 6N
N log2 N, (3.9)
where symbol “” means that asymptotically both sides have the same order. More specif-
ically, we have
C(N )
lim = a constant.
N →+∞ N log2 N
We see that the complexity of the DFT can be as low as O(N log2 N ). Since the matrix
corresponding to the inverse DFT (IDFT) is the element-wise conjugate of the matrix
corresponding to DFT, the same argument also applies to the IDFT. Hence there is an
O(N log2 N ) algorithm for the inverse discrete Fourier transform as well.
We can utilize FFT to do fast convolution, According to the following theorem.
of length N , N ∈ N:
X = {X0 , X1 , . . . , XN −1 },
Y = {Y0 , Y1 , . . . , YN −1 }.
and Y denote the discrete Fourier transform (as defined in (3.1)) of X and Y
Let X
respectively. Let X ∗ Y denote the convolution of these two sequences:
k
N −1
X ∗ Y [k] = X[i]Y [k − i] + X[i]Y [N + k − i], k = 0, 1, 2, . . . , N − 1.
i=0 i=k+1
Let X ∗ Y denote the discrete Fourier transform of the sequence X ∗ Y . We have
X Y [k],
∗ Y [k] = X[k] k = 0, 1, 2, . . . , N − 1. (3.10)
In other words, the discrete Fourier transform of a convolution of any two sequences, is equal
to the elementwise multiplication of the discrete Fourier transforms of these two sequences.
C (N ) = 2C(N ) + N + C(N )
= 3C(N ) + N
by (3.9)
N log2 N.
Hence we can do fast convolution with complexity O(N log2 N ), which is lower than order
O(N 2 ).
38 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
Definition
Coefficients of the discrete cosine transform (DCT) are equal to the inner product of a
discrete signal and a discrete cosine function. The discrete cosine function is an equally
spaced sampling of a cosine function. Note that we stay in a discrete setting. Let the
sequence
X = {X[l] : l = 0, 1, 2, . . . , N − 1}
Y = {Y [k] : k = 0, 1, 2, . . . , N − 1}.
N −1
Y [k] = Ckl X[l], 0 ≤ k ≤ N − 1,
l=0
The purpose of choosing different values for the function αi (k) is to make the DCT
an orthogonal transform, or equivalently to make the transform CN = {Ckl } k=0,1,... ,N −1 an
l=0,1,... ,N −1
orthonormal matrix.
A good reference book about DCT is Rao and Yip [121].
By choosing different values for δ1 and δ2 in (3.11), we can define four types of DCT, as
summarized in Table 3.2.
3.1. DCT AND HOMOGENEOUS COMPONENTS 39
We can view the four types of DCT as choosing the first index (which usually corresponds
to the frequency variable) and the second index (which usually corresponds to the time
variable) at different places (which can be either integer points or half integer points).
Figure 3.1 depicts this idea.
where R is a real operator: R(A) is the real part of matrix A. Symbol “•” represents an
arbitrary N × N matrix.
The other types of DCT, generally speaking, can be implemented via the 8N -point DFT.
We explain the idea for the case of type-IV DCT. For type-II and type-III DCT, a similar
idea should work. Note that in practice, there are better computational schemes.
Let us first give some mathematical notation. Let symbol IN →8N denote an insert
operator: it takes an N -D vector and expand it into an 8N -D vector, so that the subvector
at locations {2, 4, 6, . . . , 2N } of the 8N -D vector is exactly the N -D vector, and all other
40 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
II VI II VI II VI II VI II
4 I III I III I III I III I
II VI II VI II VI II VI II
second index
3 I III I III I III I III I
II VI II VI II VI II VI II
2 I III I III I III I III I
II VI II VI II VI II VI II
1 I III I III I III I III I
1 2 3 4 5
first index
Figure 3.1: Sampling of four types of DCT at the index domain. For example, both the
first index and the second index for the type-I DCT take integral values.
3.1. DCT AND HOMOGENEOUS COMPONENTS 41
elements in the 8N -D vector are zero. Equivalently, for {X[0], X[1], . . . , X[N − 1]} ∈ RN ,
Let symbol S8N →N denote the inverse operator of IN →8N . The operator S8N →N is actually
a downsampler:
S8N →N ({•, X[0], •, X[1], •, . . . , •, X[N − 1], •, . . . , •}) = {X[0], X[1], . . . , X[N − 1]},
6N
where “•” represents any complex number. If DFT8N denotes the 8N -point DFT operator
and DCTIV
N denotes the N -point type-IV DCT operator, then the following is true (note
the index of DFT starts at zero):
In order to do an N -point type-IV DCT, we can first insert the sequence into an eight-times
expanded sequence, and then apply the 8N -point DFT to the expanded vector, and then
downsample the output of DFT, and then take the real part of the downsampled result.
Note in this section, we always assume that the input is a real sequence or a real vector.
The reason that we need a much higher dimension (8 times) in DFT to implement DCT
is because DCT could be sampled at half integer points (multipliers of 1/2).
We summarize the development of fast algorithms for DCTs. Since the type-II DCT
has been used in JPEG—an international standard for image coding, compression and
transmission—research in this area has been very active. The 2-D DCT on small blocks—
for example, an 8 × 8 block or an 16 × 16 block—is of particular importance. JPEG uses
2-D DCT on an 8 × 8 block.
Let’s start with the 1-D DCT. Table 3.3 summarizes some key results. (The table is
originally from Antonio Ortega for a talk given at ISL, Stanford University.) By definition,
the 1-D N -point DCT takes N 2 multiplications and N (N − 1) additions; when N = 8,
it is equivalent to 64 multiplications and 56 additions. In 1977 [28], Chen, Smith and
42 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
Fralick developed a fast algorithm that would take N log2 N − 3N/2 + 4 multiplications and
3N/2(log2 N − 1) + 2 additions; when N = 8, their algorithm would take 16 multiplications
and 26 additions. This is a significant reduction in complexity. In 1984, Wang [139] gave an
algorithm that would take N (3/4 log2 N − 1) + 3 multiplications and N (7/4 log2 N − 2) + 3
additions; when N = 8, this is 13 multiplications and 29 additions. Also in 1984, Lee [92]
introduced an algorithm that would take N/2 log2 N multiplications and 3N/2 log2 N −N +1
additions; when N = 8, this is 12 multiplications (one less than 13) and 29 additions.
Table 3.3: Number of multiplications and additions for various fast one-D DCT/IDCT
algorithms. The number in { · } is the value when N is equal to 8.
The overall complexity of the DCT arises from two parts: multiplicative complexity
and additive complexity. The multiplicative complexity is the minimum number of nonra-
tional multiplications necessary to perform DCTs. The additive complexity, correspond-
ingly, is the minimum number of additions necessary to perform DCTs. Since it is more
complex to implement multiplications than additions, it is more important to consider mul-
tiplicative complexity. In 1987, Duhamel [62] gave a theoretical bound on the 1-D N -point
DCT: it takes at least 2N − log2 N − 2 multiplications. In 1992, Feig and Winograd [67]
extended this result to an arbitrary dimensional DCT with input sizes that are powers
of two. Their conclusion is that for L-dimensional DCTs whose sizes at coordinates are
2m1 , 2m2 , . . . , 2mL , with m1 ≤ m2 ≤ . . . mL , the multiplicative complexity is lower bounded
by 2m1 +m2 +...+mL−1 (2mL +1 − mL − 2). When L = 1, this result is the same as Duhamel’s
result [62].
In image coding, particularly in JPEG, 2-D DCTs are being used. In 1990, Duhamel
3.1. DCT AND HOMOGENEOUS COMPONENTS 43
and Guillemot [61] derived a fast algorithm for the 2-D DCT that would take 96 multipli-
cations and more than 454 additions for an 8 × 8 block. In 1992, Feig and Winograd [66]
proposed another algorithm that would take 94(< 96) multiplications and 454 additions for
an 8 × 8 block. They also mentioned that there is an algorithm that would take 86 (< 94)
multiplications, but it would take too many additions to be practical.
In practical image coding/compression, a DCT operator is usually followed by a quan-
tizer. Sometimes it is not necessary to get the exact values of the coefficients; we only need
to get the scaled coefficients. This leads to a further saving in the multiplicative complexity.
In [66], Feig and Winograd documented a scaled version of DCT that would take only 54
multiplications and 462 additions. Moreover, in their fast scaled version of the 2-D DCT,
there is no computation path that uses more than one multiplication. This makes parallel
computing feasible.
As for the DCT, the output of the 1-D discrete sine transform (DST) is the inner product
of the 1-D signal and the equally sampled sine function. The sampling frequencies are
also equally spaced in the frequency domain. Again, there are four types of DST. More
specifically, if we let
X = {X[l] : l = 0, 1, 2, . . . , N − 1}
Y = {Y [k] : k = 0, 1, 2, . . . , N − 1},
we have
N −1
Y [k] = SN (k, l)X[l], k = 0, 1, 2, . . . , N − 1,
l=0
where SN (k, l), k, l = 0, 1, 2, . . . , N − 1, are values of sine functions. For the four types of
DST, let
√1 if k = 0 or N,
bk := 2 (3.12)
1 if k = 1, . . . , N − 1;
44 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
then we have
• DST-I: SN (k, l) = 2
sin πkl
N N ;
π(k+1)(l+ 21 )
• DST-II: SN (k, l) = bk+1 N2 sin N ;
π(k+ 21 )(l+1)
• DST-III: SN (k, l) = bl+1 N2 sin N ;
π(k+ 21 )(l+ 21 )
• DST-IV: SN (k, l) = N2 sin N .
III −1
DSTII
N = (DSTN ) .
The reason that the DCT is so powerful in analyzing homogeneous signals is that it is
nearly (in some asymptotic sense) the Karhunen-Loéve transform (KLT) of some Gaussian
Markov Random Fields (GMRFs). In this section, we first describe the definition of Gaus-
sian Markov Random Fields, then argue that the covariance matrix is the key statistic for
GMRFs; for a covariance matrix, we give the necessary and sufficient conditions of diag-
onalizability of different types of DCTs; finally we conclude that under some appropriate
boundary conditions, the DCT is the KLT of GMRFs.
As we stated earlier, in this thesis, not much attention is given to mathematical rigor.
This subsubsection is organized in the following way: we start with the definition of a
random field, then introduce the definition of a Markov random field and a Gibbs random
field; the Hammersley-Clifford theorem creates an equivalence between a Markov random
field and a Gibbs random field; then we describe the definition of a Gaussian Markov random
field; eventually we argue that DCT is a KLT of a GMRF.
3.1. DCT AND HOMOGENEOUS COMPONENTS 45
Definition of a random field. We define a random field on a lattice. Let Zd denote the
d-dimensional integers, or the lattice points, in d-dimensional space, which is denoted by
Rd . The finite set D is a subset of Zd : D ⊂ Zd . For two lattice points x, y ∈ Zd , let |x − y|
denote the Euclidean distance between x and y. The set D is connected if and only if for any
x, y ∈ D, there exists a finite subset {x1 , x2 , . . . , xn } of D, n ∈ N, such that (1) |x − x1 | ≤ 1,
(2) |xi − xi+1 | ≤ 1, i = 1, 2, . . . , n − 1, and (3) |xn − y| ≤ 1. We call a connected set D a
domain. The dimension of the set D is, by definition, the number of integer points in the
set D. We denote the dimension of D by dim(D). On each lattice point in set D, a real
value is assigned. The set RD , which is equivalent to Rdim(D) , is called a state space. Follow
some conventions, we denote the state space by Ω, so we have Ω = Rdim(D) . Let F be the
σ-algebra that is generated from the Borel sets in Ω. Let P be the Lebesgue measure. The
triple (Ω, F, P) is called a random field (RF) on the domain D.
Note that we define a random field on a subset of all the lattice points.
Now we give the definition of a neighbor. Intuitively, under Euclidean distance, two
integer (lattice) points x and y are neighbors when |x − y| ≤ 1. This definition can be
extended. We define a non-negative, symmetric and translation-invariant bivariate function
N (x, y) on domain D, such that for x, y ∈ D, the function N satisfies
1. N (x, x) = 0,
2. N (x, y) ≥ 0 (non-negativity),
Any two points are called neighbors if and only if N (x, y) > 0. For example, in Euclidean
space, if we let N (x, y) = 1 when |x − y| = 1, and N (x, y) = 0 elsewhere, then we have the
ordinary definition of neighbor that is mentioned at the beginning of this paragraph.
Definition of a Markov random field. The definition of a Markov random field is based
upon conditional probability. The key idea of Markovity is that conditional probability
should depend only on neighbors. To be more precise, we need some terminology. Let ω
denote an element of Ω. We call ω a realization. Let p(ω) denote the probability density
function of ω. The p.d.f. p is associated with the Lebesgue measure P. Let ω(x) be the
value of the realization ω at the point x. For a subset A ⊂ D, suppose the values at points
46 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
in A are given by a deterministic function f (·), the conditional p.d.f. at point x ∈ D given
values in A is denoted by p{ω(x)|ω(y) = f (y), y ∈ A ⊂ D}. To simplify, we let
where the constant Z is a normalizing constant, and the function U (x, y) has the following
properties: for x, y ∈ D, we have
3. U (x, y) = 0 when points x and y are not neighbors (nearest neighbor property).
The rigorous statement and the actual proof is too long to be presented here, readers
are referred to [12, 129] for technical details.
Definition of a Gaussian Markov random field. In (3.13), if the p.d.f. p(ω) is also the
p.d.f. of a multivariate Normal distribution, then there are two consequences: first, from
the Hammersley-Clifford Theorem, it is a Markov random field; second, the random field is
also Gaussian. We call such a random field a Gaussian Markov random field (GMRF). Let
ω denote a vector corresponding to a realization ω. (The vector ω is simply a list of all the
values taken by ω.) We further suppose the vector ω is a column vector. Let Σ denote the
covariance matrix of the corresponding multivariate Normal distribution. We have
− 12 dim(D) 1 1 T −1
p(ω) = (2π) exp − ω Σ ω .
det1/2 (Σ) 2
The basic idea is that if the covariance matrix can be diagonalized by a type of DCT,
then the covariance matrix must be written as a summation of two matrices. Moreover,
the two matrices must have some particular properties. One matrix must be Toeplitz
and symmetric, or such a matrix multiplied by some diagonal matrices. The other matrix
must be Hankel and counter symmetric, or counter antisymmetric, or a matrix obtained by
shifting such a matrix and then multiplying it with a diagonal matrix. Furthermore, for an
N × N covariance matrix, the dimensionality (or the degrees of freedom) of the matrix is
no higher than N .
To be more clear, we have to use some notation. Let N denote the length of the signal,
so the size of the covariance matrix is N × N . Let bk , k = 0, 1, . . . , N − 1 be the sequence
defined in (3.12). Let B be the diagonal matrix
⎛ ⎞
b0
⎜ ⎟
⎜ b1 ⎟
⎜ ⎟
B=⎜ .. ⎟. (3.14)
⎜ . ⎟
⎝ ⎠
bN −1
Let {c0 , c1 , c2 , . . . , c2N −1 } and {c0 , c1 , c2 , . . . , c2N −1 } be two finite real-valued sequences
whose elements ci and ci satisfy a symmetric or an antisymmetric condition:
ci = c2N −i i = 0, 1, 2, . . . , 2N − 1,
ci = −c2N −i i = 0, 1, 2, . . . , 2N − 1.
Note that when i = N , the second condition implies cN = 0. Suppose C1 and C1 are the
Toeplitz and symmetric matrices
⎛ ⎞ ⎛ ⎞
c0 c1 ... cN −1 c0 c1 ... cN −1
⎜ ⎟ ⎜ ⎟
⎜ c1 c0 cN −2 ⎟ ⎜ c1 c0 cN −2 ⎟
⎜ ⎟ ⎜ ⎟
C1 = ⎜ .. .. .. ⎟, and C2 = ⎜ .. .. .. ⎟ . (3.15)
⎜ . . . ⎟ ⎜ . . . ⎟
⎝ ⎠ ⎝ ⎠
cN −1 cN −2 . . . c0 cN −1 cN −2 ... c0
3.1. DCT AND HOMOGENEOUS COMPONENTS 49
⎛ ⎞
c1 c2 ... cN −1 0
⎜ ⎟
⎜ c c3 −cN −1 ⎟
⎜ 2 ... 0 ⎟
⎜ .. .. .. .. ⎟
⎜
C2 = ⎜ . ..
. ⎟. (3.17)
. . . ⎟
⎜ ⎟
⎜ c 0 ... −c3 −c2 ⎟
⎝ N −1 ⎠
0 −cN −1 . . . −c2 −c1
Note that C2 is counter symmetric (if a{i,j} is the {i, j}th element of the matrix, then
a{i,j} = a{N +1−i,N +1−j} ,for i, j = 1, 2, . . . , N ) and C2 is counter antisymmetric (if a{i,j} is
the {i, j}th element of the matrix, then a{i,j} = −a{N +1−i,N +1−j} , for i, j = 1, 2, . . . , N ).
Let D denote a down-right-shift operator on matrix, such that
⎛ ⎞
c0 c1 ... cN −2 cN −1
⎜ ⎟
⎜ c ⎟
⎜ 1 c2 ... cN −1 cN ⎟
⎜ . . .. .. ⎟
⎜
D(C2 ) = ⎜ .. .. ..
. ⎟, (3.18)
. . ⎟
⎜ ⎟
⎜ cN −2 cN −1 . . . c4 c3 ⎟
⎝ ⎠
cN −1 cN ... c3 c2
and
⎛ ⎞
c0 c1 ... cN −2 cN −1
⎜ ⎟
⎜ c c2 cN −1 ⎟
⎜ 1 ... 0 ⎟
⎜ . . .. .. ⎟
D(C2 ) = ⎜
⎜ .
. .. ..
.
. . ⎟.
⎟ (3.19)
⎜ ⎟
⎜ c −c4 −c3 ⎟
⎝ N −2 cN −1 . . . ⎠
cN −1 0 ... −c3 −c2
Let Σ denote the covariance matrix, and let Σ1 and Σ2 denote two matrices that have
special structures that we will specify later.
We will refer to the following table (Table 3.4) in the next Theorem.
Table 3.4: Decomposition of covariance matrix that can be diagonalized by different types
of DCT.
Σ = Σ1 + Σ2 = B · C1 · B + B · D(C2 ) · B.
The matrices B, C1 , C1 , C2 , C2 , D(C2 ) and D(C2 ) are defined in equations (3.14), (3.15),
(3.16), (3.17), (3.18) and (3.19), respectively.
Moreover, the elements ci and ci are determined by the following expressions:
N −1
1 πki
ci = λk b2k cos , i = 0, 1, . . . , 2N − 1,
N N
k=0
N −1
1 π(k + 12 )i
ci = λk cos , i = 0, 1, . . . , 2N − 1,
N N
k=0
Why is this result meaningful? Because it specifies the particular structure of matrices
that can be diagonalized by DCTs. Some previous research has sought the specifications
(for example, the boundary conditions) of a GMRF so that the covariance matrix has one
of the structures that are described in Table 3.4, or has a structure close to one of them.
The research reported in [110, 111] is an example.
The 2-D DCT is an extension of the 1-D DCT. On a matrix, we simply apply a 1-D DCT
to each row, and then apply a 1-D DCT to each column. The basis function of the 2-D
DCT is the tensor product of (two) basis functions for the 1-D DCT.
The optimality of the 2-D DCT is nearly a direct extension of the optimality of the 1-D
DCT. As long as the homogeneous components of images are defined consistently with the
definition corresponding to 1-D homogeneous components, we can show that the 2-D DCT
will diagonalize the covariance matrices of the 2-D signal (images) too. Based on this, we
can generate a slogan: “The 2-D DCT is good for images with homogeneous components.”
Another way to show the optimality of the DCT is to look at its basis functions. The
properties of the basis functions represent the properties of a transform. A coefficient after a
transform is equal to the inner product of an input and a basis function. Intuitively, if basis
functions are similar to the input, then we only need a small number of basis functions to
represent the original input. Note here that the input is a signal. When the basis functions
are close to the features we want to capture in the images, we say that the corresponding
transform is good for this type of images.
To show that the 2-D DCT is good for homogeneous images, let’s look at its basis func-
tions. Figure 3.2 shows all 2-D DCT basis functions on an 8 × 8 block. From these figures,
we should say that the basis functions are pretty homogeneous, because it is impossible to
say that the statistical properties in one area are different from the statistical properties in
another.
It is interesting to review the history of the development of wavelets. The following de-
scription is basically from Deslauriers and Dubuc [60, Preface]. In 1985 in France, the first
orthogonal wavelet basis was discovered by Meyer. Shortly thereafter Mallat introduced
multiresolution analysis, which explained some of the mysteries of Meyer’s wavelet basis.
In February, 1987, in Montreal, Daubechies found an orthonormal wavelet basis that has
compact support. Next came the biorthogonal wavelets related to splines. In Paris, in May
1985, Deslauriers and Dubuc used dyadic interpolation to explain multiresolution, a term
that by then was familiar to many. Donoho and Johnstone developed the wavelet shrinkage
3.2. WAVELETS AND POINT SINGULARITIES 53
method in denoising, density estimation, and signal recovery. Wavelets and related technolo-
gies have a lot of influence in contemporary fields like signal processing, image processing,
statistic analysis, etc. Some good books that review this literature are [20, 34, 100, 107].
It is almost impossible to cover this abundant area in a short section. In the rest of the
section, we try to summarize some key points. The ones we selected are (1) multiresolution
analysis, (2) filter bank, (3) fast discrete algorithm, and (4) optimality in processing signals
that have point singularities. Each of the following subsections is devoted to one of these
subjects.
Multiresolution analysis (MRA) is a powerful tool, from which some orthogonal wavelets
can be derived.
It starts with a special multi-layer structure of square integrable functional space L2 .
Let Vj , j ∈ Z, denote some subspaces of L2 . Suppose the Vj ’s satisfy a special nesting
structure, which is:
At the same time, the following two conditions are satisfied: (1) the intersection of all
#
the Vj ’s is a null set: j∈Z Vj = ∅; (2) the union of all the Vj ’s is the entire space L2 :
j∈Z Vj = L .
2
We can further assume that the subspace V0 is spanned by functions φ(x − k), k ∈ Z,
where φ(x − k) ∈ L2 . By definition, any function in subspace V0 is a linear combination of
functions φ(x − k), k ∈ Z. The function φ is called a scaling function. The set of functions
{φ(x−k) : k ∈ Z} is called a basis of space V0 . To simplify, we only consider the orthonormal
basis, which, by definition, gives
1, k = 0,
φ(x)φ(x − k)dx = δ(k) = k ∈ Z. (3.20)
0, k = 0,
A very essential assumption in MRA is the 2-scale relationship, which is: for ∀j ∈ Z, if
function f (x) ∈ Vj , then f (2x) ∈ Vj+1 . Consequently, we can show that {2j/2 φ(2j x − k) :
k ∈ Z} is an orthonormal basis of the functional space Vj . Since φ(x) ∈ V0 ⊂ V1 , there
3.2. WAVELETS AND POINT SINGULARITIES 55
The two-scale relationship (3.21) can also be expressed in the Fourier domain:
ω ω
Φ(ω) = H( )Φ( ).
2 2
∞
$ ω
Φ(ω) = Φ(0) H( ). (3.22)
2j
j=1
Since the functional subspace Vj is a subset of the functional subspace Vj+1 , or equiva-
lently Vj ⊂ Vj+1 , we can find the orthogonal complement of the subspace Vj in space Vj+1 .
We denote the orthogonal complement by Wj , so that
Vj ⊕ Wj = Vj+1 . (3.23)
Since ψ(x) is in V1 , there must exist a real-valued sequence {g(n) : n ∈ Z} such that
ψ(x) = g(n)φ(2x − n). (3.25)
n∈Z
Actually, if relations (3.27), (3.28) and (3.29) hold, then the sequence g can only be the
reverse of the sequence h modulated by the sequence {(−1)n : n ∈ Z}. The following
theorem provides a formal description.
Theorem 3.6 (Quadratic Mirror Filter) If two sequences {h(n) : n ∈ Z} and {g(n) :
n ∈ Z} satisfy the relations listed in (3.27), (3.28) and (3.29), then the sequence {g(n) :
n ∈ Z} is uniquely determined by the sequence {h(n) : n ∈ Z} via
After we specify the sequence {h(n) : n ∈ Z}, the sequence {g(n) : n ∈ Z} is determined.
The Fourier transform of {h(n) : n ∈ Z}, which is denoted by H(ω), is determined too.
From relation (3.22), the function φ(x) is determined, so the sequence {h(n) : n ∈ Z} plays
a deterministic role in the whole design scheme. The choice of {h(n) : n ∈ Z} is determined
by the conditions that are imposed on the wavelet functions. Some such conditions are
symmetry, vanishing moments, compact support, etc. Details on how to design wavelet
functions (or, equivalently, how to choose the h sequence) are not directly relevant to this
thesis. Interested readers are referred to [34] and other books.
HIGH DOWN
PASS SAMPLER
FILTER
HIGH DOWN
SIGNAL PASS SAMPLER
FILTER
LOW DOWN HIGH DOWN
PASS SAMPLER PASS SAMPLER
FILTER FILTER
LOW DOWN
PASS SAMPLER
FILTER
LOW DOWN
PASS SAMPLER
FILTER
Figure 3.3: Illustration of a filter bank for forward orthonormal wavelet transform.
Nowadays, computers are widely used in scientific computing, and most of a signal pro-
cessing job is done by a variety of chips. All chips use digital signal processing (DSP)
technology. It is not an exaggeration to say that a technique is crippled if it does not have
a corresponding discrete algorithm. A nice feature of wavelets is that it has a fast discrete
algorithm. For length N signal, the order of computational complexity is O(N ). It can be
formulated as an orthogonal transform. Each of the wavelets can have finite support.
To explain how the discrete wavelet transform (DWT) works, we need to introduce some
notation. Consider a function f ∈ L2 . Let αkj denote a coefficient of a transform of f ; αkj
corresponds to the kth scaling function at the jth scale:
αkj = f (x)φ(2j x − k)dx.
Note that coefficients αkj and βkj are at the same scale. The following two relations can be
3.2. WAVELETS AND POINT SINGULARITIES 59
βkj = f (x)ψ(2j x − k)dx
(3.25)
= g(n) f (x)φ(2j+1 x − 2k − n)dx
n∈Z
j+1
= g(n)αn+2k . (3.32)
n∈Z
These two relations, (3.31) and (3.32), determine the two-scale relationship in the discrete
algorithm.
Vj = V0 ⊕ W0 ⊕ . . . ⊕ Wj−1 .
αj−1 = (h ∗ αj ) ↓ 2. (3.33)
60 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
Similarly, from (3.32), for any j ∈ N, the sequence β j−1 is a downsampled version of a
convolution of two sequences αj and g:
αj−1 = (h ∗ αj ) ↓ 2. (3.34)
The algorithm for the DWT arises from a combination of equations (3.33) and (3.34). We
start with αj in subspace Vj , j ∈ N. From (3.33) and (3.34) we get sequences αj−1 and β j−1
corresponding to decompositions in subspaces Vj−1 and Wj−1 , respectively. Then we further
map the sequence αj−1 into sequences αj−2 and β j−2 corresponding to decompositions
in subspaces Vj−2 and Wj−2 . We continue this process until it stops at scale 0. The
concatenation of sets α0 , β 0 , β 1 , . . . , β j−1 gives the decomposition into a union of all the
subspaces V0 , W0 , W1 , . . . , Wj−1 . The sequence α0 corresponds to a projection of function
f into subspace V0 . For l = 0, 1, . . . , j − 1, the sequence β l corresponds to a projection of
function f into subspace Wl . The set α0 ∪ β 0 ∪ β 1 ∪ . . . ∪ β j−1 is called the DWT of the
sequence αj .
The first graph in Figure 3.4 gives a scaled illustration of the DWT algorithm for se-
quence α3 . Note that we assume a sequence has finite length. The length of each block is
proportional to the length of its corresponding coefficient sequence, which is a sequence of
α’s and β’s. We can see that for a finite sequence, the length of the DWT is equal to the
length of the original sequence. The second graph in Figure 3.4 is a non-scaled depiction
for a more generic scenario.
We will explain why the wavelet transform is good at processing point singularities. The
key reason is that a single wavelet basis function is a time-localized function. In other
words, the support of the basis function is finite. Sometimes we say that wavelet basis
function has compact support. The wavelet basis is a system of functions made by dilations
and translations of a single wavelet function, together with some scaling functions at the
coarsest scale. In the discrete wavelet transform, for a signal concentrated at one point
(or, equivalently, for one point singularity) at every scale because of the time localization
property, there are only a few significant wavelet coefficients. For a length-N signal, the
number of scales is O(log N ), so for a time singularity, a DWT should give no more than
O(log N ) wavelet coefficients that have significant amplitudes. Compared with the length
3.2. WAVELETS AND POINT SINGULARITIES 61
α3
SIGNAL
@
L H @
R
@
α2 β2
@
L H @
R
@
α1 β1
@
L H@
R
@
α0 β0
j
t α JJJ
yttt %
αj−1 II β j−1
u I$
uz u
αj−2 II β j−2
I$
β j−3
z
α1 A
~}}} A
α0 β0
Figure 3.4: Illustration of the discrete algorithm for forward orthonormal wavelet transform
on a finite-length discrete signal. The upper graph is for cases having 3 layers. The width
of each block is proportional to the length of the corresponding subsequence in the discrete
signal. The bottom one is a symbolic version for general cases.
62 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
of the signal (namely N ), O(log N ) can be treated as a constant. Hence, we say that the
wavelet transform is good at processing point singularities.
Figure 3.5 illustrates a multiresolution analysis of a time singularity. The top left curve
is a function that contains a time singularity—it is generated by a Laplacian function. The
other plots in the left column show its decompositions into the coarsest scale subspace
V0 (made by dilation and translation of a scaling function) and the wavelet subspace W0
through W5 . The right column contains illustrations of the DWT coefficients corresponding
to different functional subspaces at different scales. Note that at every scale, there are only
a few significant coefficients. Compared with the length of the signal (which is N ), the total
number of the significant DWT coefficients (O(log N )) is relatively small.
signal
Wavelet Coefficients
V0
W
0
W
1
W2
W
3
W
4
W5
2-D Wavelets
The conclusion that the 2-D wavelet is good at processing point singularities in an image
is a corollary of the 1-D result.
We consider a 2-D wavelet as a tensor product of two 1-D wavelets. It is easy to observe
that the 2-D wavelet is a spatially localized function. Utilizing the idea from the filter bank
and multi-resolution analysis, we can derive a fast algorithm for 2-D wavelet transform.
The fast algorithm is also based on a two-scale relationship.
Using a similar argument as in the 1-D case, we can claim that for an N × N image,
no more than O(log N ) significant 2-D wavelet coefficients are needed to represent a point
singularity in an image. Hence the wavelet transform is good at processing point singularities
in images.
Figure 3.6 shows some 2-D wavelet basis functions. We see that wavelet functions look
like points in the image domain.
Figure 3.6: Two-dimensional wavelet basis functions. These are 32 × 32 images. The upper
left one is a tensor product of two scaling functions. The bottom right 2 by 2 images and
the (2, 2)th image are tensor products of two wavelets. The remaining images are tensor
products of a scaling function with a wavelet.
64 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
The edgelet system, defined in [50], is a finite dyadically-organized collection of line segments
in the unit square, occupying a range of dyadic locations and scales, and occurring at a range
of orientations. It is a low cardinality system. This system has a nice property: any line
segment in an image (in the system or not) can be approximated well by a few elements
from this system. More precisely, it takes at most 8 log2 (N ) edgelets to approximate any
line segment within distance 1/N + δ, where N is the size of the image and δ is a constant.
The so-called edgelet transform takes integrals along these line segments.
The edgelet system is constructed as follows:
[E1] Partition the unit square into dyadic sub-squares. The sides of subsquares are 1/2,
1/4, 1/8, . . . .
[E2] On each subsquare, put equally-spaced vertices on the boundary. The inter-distance
is prefixed and dyadic.
[E3] If a line segment connecting any two vertices is not on a boundary, then it is an edgelet.
[E4] The edgelet system is a collection of all the edgelets as defined in [E3].
Edgelets are not functions and do not make a basis; instead, they can be viewed as geometric
objects—line segments in the square. We can associate these line segments with linear
%
functionals: for a line segment e and a smooth function f (x1 , x2 ), let e[f ] = e f . Then
the edgelet transform can be defined as the mapping: f → {e[f ] : e ∈ En }, where En is the
%
collection of edgelets, and e[f ] is, as above, the linear functional e f .
Some examples of the edgelet transform on images are given in Appendix A.
3.3. EDGELETS AND LINEAR SINGULARITIES 65
The edgelet transform (defined in previous sections) is an O(N 3 log2 N ) algorithm. The
order of complexity can be reduced to O(N 2 log2 N ), if we allow a slight modification of the
original edgelet system.
Recently, a fast algorithm for calculating a discrete version of the Radon transform has
been developed. The Radon transform of a 2-D continuous function is simply a projection
of the 2-D function onto lines passing through the origin. In discrete case, people can utilize
the idea from the Fourier Slice Theorem (see Appendix B) to design a fast algorithm. It
takes the following three steps:
The so-called fast approximate edgelet transform is essentially the discrete Radon transform
of an image at various scales: we partition a unit square into subsquares having various
widths, and then apply the discrete version of the Radon transform on images limited within
each of these subsquares. This method was described in [48] and will be described in detail
in Appendix B. As we can see, the existence of this fast algorithm relies on the existence
of FFT and fractional FT.
The Riesz representers of the fast approximate edgelet transform are close to edgelets
but not exactly the same. Illustrations of some of these representers and some examples
66 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
of this transform on images are given in Appendix B, together with detailed discussion on
algorithm design issues and analysis.
In this section, we review some other widely known transforms. The first subsection is
for one-dimensional (1-D) transforms. The second subsection is for 2-D transforms. We
intend to be comprehensive, but sacrifice depth and mathematical rigor. The motivation
for writing this section is to provide us with some starting points in this field.
We briefly describe some 1-D transforms. Note that all the transforms we mention are
decomposition schemes.
As mentioned at the beginning of this chapter, there is another popular idea in designing
the 1-D transform. This idea is to map the 1-D signal into a 2-D function, which is sometimes
called a distribution. This idea is called the idea of distribution. Typically the mapping is
a bilinear transform; for example, the Wigner-Ville distribution. In this thesis, we are not
going to talk about the distribution idea. Readers can learn more about it in Flandrin’s
book [68].
We focus on the idea of time-frequency decomposition. Suppose every signal is a super-
position of some elementary atoms. Each of these elementary atoms is associated with a
2-D function on a plane, in which the horizontal axis indicates the time and the vertical axis
indicates the frequency. This two-dimensional plane is called a time frequency plane (TFP).
The 2-D function is called a time frequency representation (TFR) of the atom. Usually
the TFR of an atom is a bilinear transform of the atom. Accordingly, the same transform
of a signal is the TFR of the signal. We pick atoms such that the TFR of every atom
resides in a small cell on the TFP. We further assume that these previously mentioned cells
cover the entire TFP. The collection of these cells is called a time frequency tiling. When
a signal is decomposed as a linear combination of these atoms, the decomposition can also
be viewed as a decomposition of a function on the TFP, so we can call the decomposition
a time-frequency decomposition.
Based on Heisenberg’s uncertainty principle, these cells, which correspond to the ele-
mentary atoms, cannot be too small. There are many ways to present the Heisenberg’s
3.4. OTHER TRANSFORMS 67
uncertainty principle. One way is to say that for any signal (here a signal corresponds
to a function in L2 ), the multiplication of the variances in both the time domain and the
frequency domain are lower bounded by a constant. From a time-frequency decomposition
point of view, the previous statement means that the area of a cell for each atom can not be
smaller than a certain constant. A further result from Slepian, Pollak and Landau claims
that a signal can never have both finite time duration and finite bandwidth simultaneously.
For a careful description about them, readers are referred to [68].
1 1 1 1
Frequency
Frequency
Frequency
Frequency
0 Time 1 0 Time 1 0 Time 1 0 Time 1
(e) Chirplet (f) fan (g) Cosine Packet (h) Wavelet Packet
1 1 1 1
Frequency
Frequency
Frequency
Frequency
Figure 3.7: Idealized tiling on the time-frequency plane for (a) sampling in time domain
(Shannon), (b) Fourier transform, (c) Gabor analysis, (d) orthogonal wavelet transform, (e)
chirplet, (f) orthonormal fan basis, (g) cosine packets, and (h) wavelet packets.
The way to choose these cells determines the time-frequency decomposition scheme.
Two basic principles are generally followed: (1) the tiling on TFP does not overlap, and (2)
the union of these cells cover the entire TFP.
Time-frequency decomposition is a powerful idea, because we can idealize different trans-
forms as different methods of tiling on the TFP. Figure 3.7 gives a pictorial tour of idealized
tiling for some transforms.
68 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
Gabor Analysis
The Gabor transform was first introduced in 1946 [69]. A contemporary description of
the Gabor transform can be found in many books, for example, [68] and [100]. For a 1-D
continuous function f (t) and a continuous time variable t, the Gabor transform is
+∞
G(m, n) = f (t) · h(t − mT )e−iω0 nt dt, for integer m, n ∈ Z, (3.35)
−∞
where h is a function with finite support, also called a window function. The constants T
and ω0 are sampling rates for time and frequency. A function
is called a basis function of Gabor analysis. (Of course we assume that the set {bm,n : m, n ∈
Z} does form a basis.) Function bm,n (t) should be both time and frequency localized.
A basis function of the Gabor transform is a shifted version of one initializing function.
The shifting is operated in both time (by mT ) and frequency (by nω0 ) domains, so the
Gabor transform can be compared to partitioning the TFP into rectangles that have the
same width and height. Figure 3.7 (c) depicts an idealized tiling corresponding to Gabor
analysis.
Choices of T and ω0 determine the resolution of a partitioning. The choice of the window
function h, which appears in both (3.35) and (3.36), determines how concentrated (in both
time and frequency) a basis function bm,n is. Two theoretical results are noteworthy:
Gabor analysis is a powerful tool. It has become a foundation for the joint time frequency
analysis.
3.4. OTHER TRANSFORMS 69
Chirplets
where h is still a window function, T is the time sampling rate, and q(t) is a quadratic
polynomial of t. An example of h is the Gaussian function h(t) = e−t . Note h is a real
2
function. Suppose the general form for the quadratic polynomial is q(t) = a2 t2 + a1 t + a0 .
The instantaneous frequency is (in this case) q (t) = 2a2 t + a1 , which is linear in t.
On the TFP, a support of a chirplet is a needle-like atom. Figure 3.7(e) shows a chirplet.
An early reference on this subject is [109]. A more general but complicated way to formulate
chirplets is described in [103]. An idealistic description of the chirplet formulation in [103]
is that we consider the chirplet as a quadrilateral on the time frequency plane. Since
each corner of the quadrilateral has two degrees of freedom, the degrees of freedom under
this formulation can go up to 8. A signal associated with this quadrilateral is actually a
transform of a signal associated with a rectangle in TFP; for example, the original signal
could be a Gabor basis function. There are eight allowed transforms: (1) translation in
time, (2) dilation in time, (3) translation in frequency, also called modulation, (4) dilation
in frequency, (5) shear in time (Fourier transformation, followed by a multiplication by a
chirp, then an inverse Fourier transformation), (6) shear in frequency (multiplication by a
chirp), (7) perspective projection in the time domain, and finally (8) perspective projection
in the frequency domain.
where the constant c is a positive real number. It can be proven that a transform Wc , for
c > 0, is isometric in L2 (R). We also know that the Fourier transform is isometric in L2 (R)
(Parseval). Suppose we have an orthonormal wavelet basis W = {wj,k : j, k ∈ Z}. We
first apply a Fourier transform on the basis W , then apply an axis warping transform, then
apply the inverse Fourier transform. Because each step is an isometric transform, the result
must be another orthonormal basis in L2 (R). Thus we get a new orthonormal basis, called
a fan basis in [7]. Figure 3.7(d) gives an idealized tiling of wavelets on the time frequency
plane. Figure 3.7(f) gives an idealized tiling of the fan basis.
Cosine Packets
Here we describe cosine packets. We review some landmarks in the historic development of
time frequency localized orthonormal basis. Real-valued time frequency localized atoms not
only have some intuitive optimality, but also have been the tool to circumvent the barrier
brought by traditional Gabor analysis. The content of this subsubsection is more abundant
than its title suggests.
Collection of localized cosine and sine functions. Cosine packets are a collection of
localized cosine and sine functions. A readily observable optimality of cosine and sine
functions is that they are real-valued. Typically, an analyzed signal is a real-valued signal,
so localized cosine and sine functions are closer to a signal than these complex-valued
functions do. We also want a basis function to be localized. The localization should be in
both time and frequency domains: in both time and frequency (Fourier) domain, a function
must either have a finite support, or decays faster than any inverse of polynomials.
Orthonormal basis for L2 (R). In Gabor analysis, finding an orthonormal basis for L2 (R)
was a topic that has been intensively studied. The Balian-Low theorem says that if we choose
a time-frequency atom following (3.36), then the atom cannot be simultaneously localized
in both time domain and frequency domain—the atom must have infinite variance either in
time or in frequency. To circumvent this barrier, Wilson in 1987 [34, page 120] proposed a
basis function that has two peaks in frequency. In [34], Daubechies gives a construction of
an orthonormal basis for L2 (R). In her construction, basis functions have exponential decay
in both time and frequency. From [34], the key idea to ideal time-frequency localization and
orthonormality in the windowed Fourier framework is to use sine and cosine rather than
complex exponentials. This is a not-so-obvious optimality of using sine and cosine functions
in basis.
3.4. OTHER TRANSFORMS 71
Folding. Another way to obtain localized time-frequency basis, instead of using the
constructive method in Daubechies [34], is to apply folding to an existing orthonormal
basis for periodic functions on a fixed interval. The idea of folding is described in [141]
and [142]. The folding operation is particularly important in discrete algorithms because if
the basis function is from the folding of a known orthonormal basis function (e.g., Fourier
basis function), since the coefficient of the transform associated with the basis is simply the
inner product of the signal with the basis function, we can apply the adjoint of the folding
operator to the signal, then calculate its inner product with the original basis function. If
there is a fast algorithm to do the original transform (for example, for DFT, there is a FFT),
then there is a fast algorithm to do the transform associated with the new time-frequency
atoms.
Best orthogonal basis. The best orthogonal basis (BOB) algorithm is proposed to achieve
the objective just mentioned. The key idea is to assign a probability distribution on a
class of signals. Usually, the class of signals is a subspace of L2 (R). Without loss of
generality, we consider signals that are restricted on the interval [0, 1). Suppose we consider
only functions whose support are dyadic intervals that have forms [2−l (k − 1), 2−l k), for
l ∈ N, k = 1, 2, . . . , 2l . At each level l, for the basis functions whose support are intervals
like [2−l (k − 1), 2−l k), the associated coefficients can be calculated. We can compute the
entropy of these coefficients. For an interval [2−l (k−1), 2−l k), we can consider its immediate
two subintervals: [2−l−1 (2k − 2), 2−l−1 (2k − 1)) and [2−l−1 (2k − 1), 2−l−1 2k). (Note [2−l (k −
1), 2−l k) = [2−l−1 (2k−2), 2−l−1 (2k−1))∪[2−l−1 (2k−1), 2−l−1 2k).) Note in the discrete case,
a set of the coefficients associated with an interval [2−l (k − 1), 2−l k) and a set of coefficients
associated with intervals [2−l−1 (2k −2), 2−l−1 (2k −1)) and [2−l−1 (2k −1), 2−l−1 2k) have the
same cardinality. In fact, there is a binary tree structure. Each dyadic interval corresponds
to a node in the tree. Suppose a node in the tree corresponds to an interval I. Two
subsequent nodes in the tree should correspond to the two subintervals of the interval I.
72 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
The idea of the BOB is the following: if the entropy associated with the two subintervals is
lower than the entropy associated with the larger dyadic interval, then we choose to divide
the larger dyadic interval. This process is repeated until there is no more partitioning to
do. The final partition of the unit interval [0, 1) determines a subset of localized cosine
functions. It’s not hard to observe that this subset of local cosine functions makes an
orthonormal basis. For a given signal, instead of computing the entropy of a statistical
distribution, we can compute the empirical entropy, which is just the entropy of a set of
coefficients. Note that the coefficients are from a transform of the signal. Based on this
empirical entropy, we can use the BOB algorithm to choose a subset of coefficients, and
correspondingly a subset of basis functions that form a basis.
Complexity. Suppose the length of the signal is N , for N ∈ N. Note that when we take
cosine functions as an orthonormal basis over an interval, the “folding” transform gives
another orthonormal system whose elements are functions that are defined on the entire
real axis, but are localized (have finite support). For every dyadic interval, we have such
an orthonormal system. Combining all these systems, we actually get the cosine packets
(CP). Suppose dyadic intervals having the same length is a cover of the entire real axis. We
say that CP elements associated with these intervals are at the same scale. For length-N
signals, we have log2 N scales. To calculate the CP coefficients at one scale, we can first
apply “folding” to the data, then apply a DCT. The folding is an O(N ) operation. The
discrete cosine transform has O(N log N ) complexity. Since we have log2 N scales, the total
complexity of computing CP coefficients is O(N log2 N ).
Wavelet Packets
In MRA, instead of dividing the low-frequency part, or a scaling space (a space spanned
by scaling functions), we can apply a quadratic mirror filter to divide the high-frequency
part, or a wavelet space (which is spanned by the wavelet functions). This idea unveils
developments of wavelet packets. Equivalently, in the filter banks algorithm, at each step,
instead of deploying a pair of filters after an LPF, we can deploy a pair of filters after an
HPF. By doing this, we partition the high-frequency part. Since a wavelet space contains
high-frequency components, partition in a high-frequency part is equivalent to partition in
a wavelet space.
In every step, we can adaptively choose to deploy a pair of filters either after an LPF
or an HPF. In this framework, all possible sets of coefficients have a binary tree structure.
3.4. OTHER TRANSFORMS 73
We can apply the idea of BOB, which is described in the previous subsubsection, to select
an orthonormal basis.
Computational complexity. Computation of the coefficients at every step involves merely
filtering. The complexity of filtering a length-N signal with a finite-length filter is no more
than O(N ), so in each step, the complexity is O(N ). Since the wavelet packets can go
from scale 1 down to scale log2 N . The total complexity of calculating all possible wavelet
packets coefficients is O(N log2 N ).
To capture some key features in an image, different transforms have been developed. The
main idea of these efforts is to construct a basis, or frame, such that the basis functions, or
the frame elements, have the interesting features. Some examples of the interesting features
are anisotropy, directional sensitivity, etc.
In the remainder of this subsection, we introduce brushlets and ridgelets.
Brushlets
The brushlet is described in [105]. A key idea is to construct a windowized smooth orthonor-
mal basis in the frequency domain; its correspondent in the time domain is called brushlets.
By constructing a “perfectly” localized basis in the frequency domain, if an original signal
has a peak in the frequency domain, then the constructed basis, brushlets, tends to capture
this feature. A 2-D brushlet is a tensor product of two 1-D brushlets. A nice property of a
2-D brushlet is that it is an anisotropic, directionally sensitive, and spatially localized basis.
More details are in Meyers’ original paper [105].
Ridgelets
f (u · x − b)
74 CHAPTER 3. IMAGE TRANSFORMS AND IMAGE FEATURES
is a ridge function, where u denotes a constant vector, b denotes a scalar constant, and
u · x is the inner product of vector u and vector x. This idea is originally from neural
networks. When a function f is a sigmoid function, the function f (u · x − b) is a single
neuron in neural networks. For Candès’ ridgelets, a key idea is to cleverly choose some
conditions for f ; for example, an admissibility condition, so that when a function f satisfies
these conditions, we not only are able to have a continuous representation based on ridge
functions that have the form f (a linear term), but also can construct a frame for a compact-
supported square-integrable functional space. The latter is based on carefully choosing a
spatial discretization on the (u, b) plane. Note that a compact-supported square-integrable
functional space should be a subspace of L2 (R). Candès [21] proves that his ridgelet system
is optimal in approximating a class of functions that are merely superpositions of linear
singularities. A function that is a superposition of linear singularities is a function that is
smooth everywhere except on a few lines. A 2-D example of this kind of function is a half
dome: for x = (x1 , x2 ) ∈ R2 , f (x) = 1{x1 >0} e−x1 −x2 . More detailed discussion is in [21].
2 2
3.5 Discussion
In the design of different transforms, based on what we have seen, the following three
principles are usually followed:
give only a few coefficients, so the result should be sparse. When this happens, we
say that the transform T has captured features in the signal f . Hence the transform
T is ideal at processing the signal f .
In this chapter, we have presented graphs of different basis functions for different
transforms. We intend to show what kind of features these transforms may capture.
• Fast algorithm. Since most of contemporary signal processing is done in digital format,
a continuous transform is unlikely to make a strong impact unless it has a fast discrete
algorithm. Some examples are the continuous Fourier transform and the continuous
wavelet transform. Both of them have fast discrete correspondents. A fast algorithm
usually means that the complexity should not be higher than a product of the original
size (which typically is N for a 1-D signal and N 2 for a 2-D N by N image) and log N
(or a polynomial of log N ).
3.6 Conclusion
We have reviewed some transforms. A key point is that none of these transform is ideal
for processing images with multiple features. In reality, an image tends to have multiple
features.
Two consequential questions are whether and how we can combine multiple transforms,
so that we can simultaneously take advantage of all of them. The remainder of this thesis
will try to answer these questions.
Note that since images are typically very large, efficient numerical algorithms are crucial
for determining if a scheme is successful or not. Readers may notice that the remainder of
this thesis is very computationally oriented.
3.7 Proofs
Suppose amn is the mn-th component of the matrix C. The proof is based on the following
fact:
N −1
N −1
1 (C is a circulant matrix) 1
amk √ e−i N kl ck √ e−i N (m+k)l
2π 2π
=
k=0
N k=0
N
N −1
1 −i 2π ml
ck e−i N kl ,
2π
= √ e N
N k=0
where 0 ≤ m, l ≤ N − 1 and i2 = −1. Based on this, one can easily verify that the Fourier
series {e−i N kl : k = 0, 1, 2, . . . , N − 1}, for l = 0, 1, 2, . . . , N − 1, are the eigenvectors of the
2π
√ √ √
matrix C and the Fourier transform of the sequence { N c0 , N c1 , . . . , N cN −1 } are the
eigenvalues.
Suppose the covariance matrix Σ can be diagonalized by a type of DCT. If the corresponding
• , then we have C • Σ(C • )T = Ω , where Ω
DCT matrix is CN N N N N is a diagonal matrix
• )T Ω C • . From (3.11),
ΩN = diag{λ0 , λ1 , . . . , λN −1 }. Equivalently, we have Σ = (CN N N
3.7. PROOFS 77
N −1
2 π(k + δ1 )(i + δ2 ) π(k + δ1 )(j + δ2 )
Σij = λk α12 (k)α2 (i)α2 (j)
cos cos
N N N
k=0
N −1 & '
1 2 π(k + δ1 )(i − j) π(k + δ1 )(i + j + 2δ2 )
= α2 (i)α2 (j) λk α1 (k) cos + cos .
N N N
k=0
1 N −1 2 π(k+δ1 )(i−j)
In the above equation, let the first term N k=0 λk α1 (k) cos N be the ijth
1 N −1 2 π(k+δ 1 )(i+j+2δ2)
element of matrix Σ1 , and let the second term N k=0 λk α1 (k) cos N be the
ijth element of matrix Σ2 . It is easy to see that the matrix Σ1 is Toeplitz, and the matrix
Σ2 is Hankel. Inserting different values of δ1 and δ2 , we establish the Theorem 3.5.
From (3.27), note that all the even terms in H(z)H(z −1 ) have zero coefficient except the
constant term. Hence we have
From (3.29), all the even terms in H(z)G(z −1 ) have zero coefficients, so we have
Equations (3.38) (3.39) and (3.40) together are equivalent to the matrix multiplication
H(z) H(−z) H(z −1 ) H(−z −1 ) 2 0
= .
G(z) G(−z) G(z −1 ) G(−z −1 ) 0 2
If both sequence {h(n) : n ∈ Z} and sequence {g(n) : n ∈ Z} have finite length (thinking of
the impulse responses of FIR filters), then the polynomial H(z)G(−z) − G(z)H(−z) must
have only one nonzero term. Since if the polynomial H(z)G(−z) − G(z)H(−z) has more
than two terms and simultaneously has finite length, (3.41) can never be equal to a constant.
We know that H(z)G(−z) − G(z)H(−z) only has one nonzero term. The polynomial
H(z −1 ) + H(−z −1 ) can only have terms with even exponents and polynomial G(−z) − G(z)
can only have terms with odd exponents, so from the previous equation, there must exist
an integer k, such that
and
Similarly, we have
Consequently,
This chapter is about combined image representation and sparse decomposition. First, in
Section 4.1, we discuss the motivation for using combined image representation, the key be-
ing that we can obtain benefits from different image transforms. Because we have combined
representations, we have an overcomplete system, or an overcomplete dictionary. Section
4.2 surveys research developments in finding sparse decompositions in an overcomplete dic-
tionary. Section 4.3 explains the optimality of using the minimum 1 norm decomposition
and explains the formulation we used in this thesis. Section 4.4 explains how we use La-
grange multipliers to transform a constrained optimization problem into an unconstrained
optimization problem. Section 4.5 is about how to choose the parameters in our method.
Section 4.6 points out that a homotopic method converges to the minimum 1 norm decom-
position. Section 4.7 describes the Newton method. Section 4.8 surveys existing methods
and softwares and explains some advantages of our approach. Section 4.9 gives a preview
about the importance of iterative methods and why we use them. Section 4.10 gives more
detail to the numerical solution of the problem. Finally, in Section 4.11, we make some
general remarks. Section 4.12 contains all relevant proofs in this chapter.
Recently, many new methods for signal/image representation have been proposed, including
wavelets, wavelet packets, cosine packets, brushlets, edgelets, and ridgelets. Typically, each
of these is good for a specific class of features, but not for others. For example, for 1-D
signal representation, wavelets are effective at representing signals made by impulses—where
81
82 CHAPTER 4. COMBINED IMAGE REPRESENTATION
“effective” means that there are only a few coefficients that have large amplitudes—while
not effective at representing oscillatory signals. At the same time, the Fourier basis is good
at representing oscillatory signals, but not at representing impulses. For 2-D images, 2-D
wavelets are effective at representing point singularities and patches; edgelets are effective
at representing linear singularities [50]. Different transforms are effective at representing
different image features. An image is usually made of several features. Combining several
transforms, we have more flexibility, hopefully enabling sparse representation.
y = T1 x1 + T2 x2 , (4.1)
where T1 and T2 are matrices whose columns are vectorized basis functions of the 2-D DCT
and the 2-D wavelet transform, respectively. x1 and x2 are coefficient vectors. y, x1 , x2 ∈
2
RN . If there is only T1 or T2 , it takes an O(N 2 ) or O(N 2 log N ) algorithm to get the
coefficient x1 or x2 . But if we want to find a sparse solution to the overcomplete system (4.1),
say, find the minimum 1 norm solution to (4.1), then we need to solve linear programming
(LP) problem,
− −
minimize eT (x+ +
1 + x1 + x2 + x2 ),
− + −
1 , x1 , x2 , x2 ≥ 0,
subject to x+
y = T1 (x+ − x− ) + T2 (x+ −
2 − x2 ),
⎛ 1⎞ 1
1
⎜ ⎟
⎜ 1 ⎟
⎜ ⎟
e = ⎜ ⎟,
⎜ 1 ⎟
⎝ ⎠
1
which becomes much more complicated. We can solve it with the simplex method or the
4.2. SPARSE DECOMPOSITION 83
interior point method. As we know, solving this LP problem in general takes more time
than doing a 2-D DCT or a 2-D wavelet transform.
Additional information about why we choose a minimum 1 norm solution is in Section
4.2 and Section 4.3.
There is by now extensive research on finding sparse decompositions. These methods can
be roughly classified into three categories:
1. greedy algorithms,
3. special structures.
with some recent advances in theoretical study, show that a global optimization algorithm
like BP is more stable in recovering the original sparse decomposition, if it exists. But BP
is a computationally intensive method. The remainder of this thesis is mainly devoted to
overcoming this barrier. In the next section, Section 4.3, we first explain the optimality of
the minimum 1 norm decomposition and then give our formulation, which is a variation
of the exact minimum 1 norm decomposition. This formulation determines the numerical
problem that we try to solve.
The minimum 1 norm solution means that in a decomposition y = Φx, where y is the desired
signal/image, Φ is a flat matrix with each column being an atom from an overcomplete
dictionary and x is a coefficient vector, we pick the one that has the minimum 1 norm
(x1 ) of the coefficient vector x.
A heuristic argument about the optimality of the minimum 1 norm decomposition is
that it is the “best” convexification of the minimum 0 norm problem. Why? Suppose we
consider all the convex functions that are supported in [−1, 1] and upper bounded by the
0 norm function and we solve
⎧
⎪
⎨ f (x) ≤ x0 ,
⎪
maximize f (x), subject to x∞ ≤ 1,
x ⎪
⎪
⎩ f is convex.
Another recent advance in theory [43] has given an interesting result. Suppose that the
overcomplete dictionary we consider is a combination of two complementary orthonormal
bases. By “complementary” we mean that the maximum absolute value of the inner product
of any two elements (one from each basis) is upper bounded by a small value. Suppose the
observed signal y is made by a small number of atoms in the dictionary. Solving the
minimum 1 norm problem will give us the same solution as solving the minimum 0 norm
problem. More specifically, if a dictionary is made by two bases—Dirac basis and Fourier
√
basis—and if the observation y is made by fewer than N /2 atoms from the dictionary,
where N is the size of the signal, then the minimum 1 norm decomposition is the same as
the minimum 0 norm decomposition.
x0 denotes the number of nonzero elements in the vector x. The 0 norm is generally
regarded as the measure of sparsity. So the minimum 0 norm decomposition is generally
regarded as the sparsest decomposition. Note the minimum 0 norm problem is a combina-
torial problem and in general is NP hard. But the minimum 1 norm problem is a convex
optimization problem and can be solved by some polynomial-time optimization methods,
86 CHAPTER 4. COMBINED IMAGE REPRESENTATION
for example, linear programming. The previous result shows that we can attack a com-
binatorial optimization problem by solving a convex optimization problem; if the solution
satisfies certain conditions, then the solution of the convex optimization problem is the same
as the solution of the combinatorial optimization problem. This gives a new possibility to
solve an NP hard problem.
More examples of identical minimum 1 norm decomposition and minimum 0 norm
decomposition are given in [43].
An exact minimum 1 norm problem is
We explain how to solve (P1 ) based on some insights from the interior point method. One
key idea is to select a barrier function and then minimize the sum of the objective function
and a multiplication of a positive constant and the barrier function. When the solution x
approaches the boundary of the feasible set, the barrier function becomes infinite, thereby
guaranteeing that the solution is always within the feasible set. Note that the subsequent
optimization problem has become a nonconstrained optimization problem. Hence we can
apply some standard methods—for example, the Newton method—to solve it.
A typical interior point method uses a logarithmic barrier function [113]. The algorithm
in [25] is equivalent to using an 2 penalty function. Since the feasible set in (P1 ) is the whole
Euclidean space, the demand of restricting the solution in a feasible set is not essential. We
actually solve the following problem
N
where λ is a scalar parameter and ρ is a convex separable function: ρ(x) = i=1 ρ̄(xi ), where
4.4. LAGRANGE MULTIPLIERS 87
ρ( ^c)
Feasible Area
1_
−
λ
2
||y||
Distortion
2
||y- Σ c i φ i ||
i=1,...,N
ρ̄ is a convex 1-D function. Note this is the idea of quadratic penalty-function method in
solving the following optimization problem with exact constraints:
To better explain the connections we raise here, we introduce a concept called Quasi-
sparsity & distortion curve. “Quasi-sparsity” refers to the small values in the quasi-sparsity
measurement ρ(x). Figure 4.1 gives a depiction. The horizontal axis is the distortion
measure y − Φx22 . The vertical axis is a measure of quasi-sparsity, in our case ρ(x).
Allowing some abuse of the terminology “sparsity”, we call this plane a distortion-sparsity
(D-S) plane. If there exists x such that (u, v) = (y − Φx22 , ρ(x)), then we say the point
(u, v) on the D-S plane is feasible. We know y − Φx22 is a quadratic form. When ρ(x) is a
convex function of x, all the feasible points on the D-S plane form a convex set; we call this
convex set a feasible area. For a fixed value of distortion, there is a minimum achievable
value for ρ(x). If all the points like this form a continuous curve, we call it a Quasi-Sparsity
88 CHAPTER 4. COMBINED IMAGE REPRESENTATION
& Distortion (QSD) curve. Actually, the QSD curve is the lower boundary of the feasible
area, as shown in Figure 4.1. Note in the figure, the notation i=1,... ,N ci φi serves the same
meaning as Φx in the above text.
A noteworthy phenomenon is that for fixed λ, the corresponding (u, v) point given by
the solution of (4.2) is actually the tangent point of a straight line having slope −1/λ with
the QSD curve. The tangent is the leftmost straight line having slope −1/λ and intersecting
with the feasible area. Moreover, the QSD curve is the pointwise upper bound of all these
tangents. We can see a similar argument on the rate and distortion curve (R&D curve) in
Information Theory [10, 33].
There are two limiting cases:
1. When the distortion y − Φx22 is zero, the QSD curve intersects with the vertical
axis. The intersection point, which has the coordinates (0, ρ(ĉ)), is associated with
the solution ĉ to the exact constraint problem as in (4.3).
2. When the measure of sparsity ρ(x) is zero, because ρ(x) is convex, nonnegative, and
symmetric about zero, we may think of x as an all-zero vector. Hence the distortion
is equal to y2 . The corresponding point on the D-S plane is (y2 , 0), and it is the
intersection of the QSD curve with the horizontal axis.
1 −γ|x| 1
ρ̄(x, γ) = |x| + e − , for γ > 0,
γ γ
where γ is a controlling parameter. Note when γ → +∞, ρ̄(x, γ) → |x|. The derivatives of
ρ̄ have the form:
∂ 1 − e−γx , x ≥ 0,
ρ̄(x, γ) = (4.4)
∂x −1 + eγx , x ≤ 0,
4.5. HOW TO CHOOSE ρ AND λ 89
2.5 1
2 0.5
ρ’(x)
ρ(x)
1.5 0
1 −0.5
0.5 −1
0 −1.5
−3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3
x x
4
ρ’’(x)
0
−3 −2 −1 0 1 2 3
x
Figure 4.2: Function ρ̄ (as ρ in figures) and its first and second derivatives. ρ̄, ρ̄ and ρ̄
are solid curves. The dashed curve in figure (a) is the absolute value. The dashed curve in
figure (b) is the signum function.
and
∂2
ρ̄(x, γ) = γe−γ|x| . (4.5)
∂2x
λ
ΦT (y − Φ
x)∞ ≤ . (4.6)
2
90 CHAPTER 4. COMBINED IMAGE REPRESENTATION
0.2
0.18
0.16
160
80
0.14
40
0.12
20
0.1
0.08 10
0.06
0.04
0.02
0
−0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2
A way to interpret the above result is that for the residual r = y − Φx̂, the maximum
amplitude of the analysis transform of the residual ΦT r∞ is upper bounded by λ2 . Hence
if λ is small enough and if ΦT is norm preserving—each column of Φ has almost the same l2
, Φ
norm—then the deviation of the reconstruction based on x x, from the image y is upper
λ
bounded by a small quantity (literally 2) at each direction given by the columns of Φ. So
if we choose a small λ, the corresponding reconstruction cannot be much different from the
original image.
After choosing ρ and λ, we have the Hessian and gradient of the objective function (later
denoted by f (x)). For fixed vector xk = (xk1 , xk2 , . . . , xkN )T ∈ RN , the gradient at xk is
⎛ ⎞
ρ̄ (xk1 )
⎜ .. ⎟
g(xk ) = −2ΦT y + 2ΦT Φxk + λ ⎜
⎝ . ⎟,
⎠ (4.7)
ρ̄ (xN ) k
where xki is the i-th element of vector xk , i = 1, 2, . . . , N , and the Hessian is the matrix
⎛ ⎞
ρ̄ (xk1 )
⎜ ⎟
H(xk ) = 2ΦT Φ + λ ⎜
⎝
..
. ⎟,
⎠ (4.8)
ρ̄ (xkN )
4.6 Homotopy
Based on the previous choice of ρ and ρ̄, when the parameter γ goes to +∞, ρ(x, γ) =
N
i=1 ρ̄(xi , γ) goes to function x1 . Note our ultimate goal is to solve the minimum
1
norm problem (P1 ). Considering the Lagrangian multiplier method, for a fixed λ, we solve
When γ goes to +∞, will the solution of the problem in (4.2) converge to the solution of
the problem in (4.9)?
We prove the convergence under the following two assumptions. The first one is easy to
satisfy. The second one seems too rigorous and may not be true for many situations that we
92 CHAPTER 4. COMBINED IMAGE REPRESENTATION
are interested in. We suspect that the convergence is still true when the second assumption
is relaxed. We leave the analysis for future research.
Assumption 1 For given y, Φ and λ, the solution to the problem in (4.9) exists and is
unique.
Since the objective function is convex and its Hessian, as in (4.8), is always positive definite,
it is easy to prove that the above assumption is true in most cases.
Theorem 4.1 For fixed γ, let x(γ) denote the solution to problem (4.2). Let x denote the
solution to problem (4.9). If the previous two assumptions are true, we have
From the above result, we can apply the following method. Starting with a small γ1 , we
get the solution x(γ1 ) of problem (4.2). Next we choose γ2 > γ1 , set x(γ1 ) as the initial guess,
and apply an iterative method to find the solution (denoted by x(γ2 )) of problem (4.2). We
repeat this process, obtaining a parameter sequence γ1 , γ2 , γ3 , . . . , and a solution sequence
x(γ1 ), x(γ2 ), x(γ3 ), . . . . From Theorem 4.1, the sequence {x(γi ), i = 1, 2, . . . } converges to
the solution of problem (4.9). This method may save computing time because when γ is
small, an iterative method takes a small number of steps to converge. After finding an
approximate solution by using a small valued γ, we then take it as a starting point for an
iterative method to find a more precise solution. Some iterations may be saved in the early
stage.
4.7. NEWTON DIRECTION 93
where βi is a damping parameter (βi is chosen by line search to make sure that the value of
the objective function is reduced), n(x(i) ) is the Newton direction, a function of the current
guess x(i) . Let f (x) = y − Φx22 + λρ(x) denote the objective function in problem (4.2).
The gradient and Hessian of f (x) are defined in (4.8) and (4.7). The Newton direction at
x(i) satisfies
( )
H(x(i) ) · n(x(i) ) = −g(x(i) ). (4.12)
This is a system of linear equations. We choose iterative methods to solve it, as discussed
in the next section and also the next chapter.
1 1
minimize cT xo + γxo 2 + p2 subject to Axo + δp = y, xo ≥ 0,
xo
2 2
where
column vector;
• xo = (xT+ , xT− )T , where x+ and x− are the positive and negative part of the column
vector x that is the same as the “x” in (4.2): x = x+ − x− , x+ ≥ 0 and x− ≥ 0;
• A = [Φ, −Φ] and Φ is the same matrix as the one specified in (4.2);
The main effort of this approach is to solve the following system of linear equations [27,
equation (6.4), page 56]
where r, v and t are column vectors given in CDS paper, D = (X −1 Z + γ 2 I)−1 , X and
Z are diagonal matrices composed from primal variable xo and dual slack variable z. (For
more specific description, we refer readers to the original paper.)
Recall in our approach, to obtain the Newton direction, we solve
* T +
2Φ Φ + λρ (xk ) n = 2ΦT (y − Φxk ) − λρ (xk ), (4.14)
where ρ (xk ) = diag{ρ̄ (xk1 ), . . . , ρ̄ (xkN )}, ρ (xk ) = (ρ̄ (xk1 ), . . . , ρ̄ (xkN ))T and n is the
desired Newton direction that is equivalent to n(x(k) ) in (4.12).
After careful examination of the matrices on the left hand sides of both (4.13) and (4.14),
we observe that both of them are positive definite matrices having the same size. By this
we say that the two approaches are similar. But it is hard to say which one may have a
better eigenvalue distribution than the other.
When they are close to the solution, both algorithms become slow to converge. Note
that the system of linear equations in (4.13) is for a perturbed LP problem. If both γ and
δ are extremely small, the eigenvalue distribution of the matrix ADAT + δ 2 I can be very
bad, so that iterative methods converge slowly. The same discussion is also true for the
matrix on the left hand side of (4.14). When the vector xk has nonzero and spread entries,
the eigenvalue distribution of the positive definite matrix 2ΦT Φ + λρ (xk ) could be bad for
any iterative methods.
The appealing properties of our approach are:
4.9. ITERATIVE METHODS 95
Chen, Donoho and Saunders [26] have successfully developed software to carry out
numerical experiments with problems having the size of thousands by tens of thousands.
We choose to implement a new approach in the hope that it simplifies the algorithm and
hopefully increase the numerical efficiency. A careful comparison of our new approach with
the one previously implemented will be an interesting future research topic. At the same
time, we choose LSQR instead of CG. (Chen et al. showed how to use LSQR when δ > 0
in (4.13), but their experiments used CG.)
4.11 Discussion
In statistics, we can find the same method being used in model selection, where we choose
a subset of variables so that the model is still sufficient for prediction and inference. To be
more specific, in linear regression models, we consider
y = Xβ + ε,
where y is the response, X is the model matrix with every column being values of a variable
(predictor), β is the coefficient vector, and ε is a vector of IID random variables. Model
selection in this setting means choosing a subset of columns of X, X (0) , so that for most
of the possible responses y, we have y ≈ X (0) β (0) , where β (0) is a subvector of β with
locations corresponding to the selected columns in X (0) . The difference (or prediction
error), y − X (0) β (0) , is negligible in the sense that it can be interpreted as a realization of
the random noise vector ε.
Typically, people use penalized regression to select the model. Basically, we solve
which is exactly the problem we encountered in (4.2). After solving problem (P R), we can
pick the ith column in X if βi has a significant amplitude. When ρ(β) = β1 , the method
(P R) is called LASSO by R. Tibshirani [133] and Basis Pursuit by Chen et al [27]. When
4.11. DISCUSSION 97
ρ(β) = β22 , the method (P R) is called ridge regression by Hoerl and Kennard [81, 80].
An ideal measure of sparsity is usually nonconvex. For example, in (4.2), the number of
nonzero elements in x is the most intuitive measure of sparsity. The 0 norm of x, x0 , is
equal to the number of nonzero elements, but it is not a convex function. Another choice
of measure of sparsity is the logarithmic function; for x = (x1 , . . . , xN )T ∈ RN , we can
have ρ(x) = N i=1 log |xi |. In sparse image component analysis, another nonconvex sparsity
measure is used: ρ(x) = N 2
i=1 log(1 + xi ) [53].
Generally speaking, a nonconvex optimization problem is a combinatorial optimization
problem, and hence it is NP hard. Some discussion about how to use reweighting methods
to solve a nonconvex optimization problem is given in the next subsection.
Sometimes, a reweighted iterative method can be used to find a local minimum for a non-
convex optimization problem. Let’s consider the following problem:
N
(LO) minimize log |xi |, subject to y = Φx;
x
i=1
N
(LOλ ) minimize y − Φx22 +λ log(|xi | + δ).
x
i=1
N
|xi |
(RIA) x(k+1) = argmin (k)
, subject to y = Φx;
i=1 |xi | + δ
x
1
More precisely, (LOλ ) is the Lagrangian multiplier version of the following optimization problem:
N
minimize log(|xi | + δ), subject to y − Φx ≤ ε.
x
i=1
N
|xi |
(RIAλ ) (k+1)
x = argmin y − Φx22 +λ (k)
.
i=1 |xi | + δ
x
(k+1) (k)
|xi − xi | → 0, as k → +∞.
(k)
Theorem 4.3 If the sequence {xi , k = 1, 2, 3, . . . } generated by (RIAλ ) converges, it
converges to a local minimum of (LOλ ).
Some related works can be found in [36, 106]. There is also some ongoing research, for
example, the work being carried out by Boyd, Lobo and Fazel in the Information Systems
Laboratories, Stanford.
4.12 Proofs
1
λ∇ρ(x) = ΦT (y − Φx) . (4.15)
2
N
Recall ρ(x) = i=1 ρ̄(xi ). It’s easy to verify that ∀xi , |ρ̄ (xi )| ≤ 1. Hence for the left-hand
side of (4.15), we have 12 λ∇ρ(x)∞ ≤ λ2 . The bound (4.6) follows. 2
4.12. PROOFS 99
We will prove convergence first, and then prove that the limiting distribution is x .
When γ takes all the real positive values, (x(γ), γ) forms a continuous and differentiable
trajectory in RN +1 . Let xiγ denote the ith element of x(γ). By the first-order condition, we
have
⎛ ⎞
∂ ρ̄(x1γ ,γ)
⎜ ∂x ⎟
⎜ .. ⎟
0 = −2ΦT (y − Φx(γ)) + λ ⎜ . ⎟.
⎝ ⎠
∂ ρ̄(xN
γ ,γ)
∂x
d
Taking dγ on both sides, we have
⎛ ⎞
∂ 2 ρ̄(x1γ ,γ) dx1γ ∂ 2 ρ̄(x1γ ,γ)
+
⎜ ∂γ∂x dγ ∂x∂x ⎟
dx(γ)
T ⎜ .. ⎟
0 = 2Φ Φ + λ⎜ . ⎟.
dγ ⎝ ⎠
∂ 2 ρ̄(xN
γ ,γ) dxN 2 N
γ ∂ ρ̄(xγ ,γ)
∂γ∂x + dγ ∂x∂x
Since
∂ 2 ρ̄(x, γ)
= xe−γ|x| ,
∂γ∂x
and
∂ 2 ρ̄(x, γ)
= γe−γ|x| ,
∂x∂x
we have
⎛ ⎞
dx1γ
x1γ e−γ|xγ | −γ|x1γ |
1
+ dγ γe
⎜ ⎟
dx(γ) ⎜ .. ⎟
−2Φ ΦT
= λ⎜ . ⎟.
dγ ⎝ ⎠
−γ|xγ | dxN −γ|xN
γ |
N
xN
γ e + γ
dγ γe
dx(γ)
Based on Assumption 2, suppose that the kth entry of vector dγ is negative and takes
the maximum amplitude of the vector:
, ,
, dx(γ) , dxk
, , = − γ.
, dγ , dγ
∞
100 CHAPTER 4. COMBINED IMAGE REPRESENTATION
Hence
, ,
xkγ e−γ|xγ |
k
, dx(γ) , dxk
, , =− γ ≤
, dγ ,
+ γe−γ|xγ |
k
dγ T
∞ 1+ 2(Φ Φ)kk
1 xkγ −γ|xk |/2
≤ √ . √ e γ
T
2 2 1+ (Φ Φ)kk γ
1 2 −1
≤ √ . 3/2
e .
T
2 2 1+ (Φ Φ)kk γ
The integration of the last term in the right-hand side of the above inequality is finite, so
the integration of dx(γ)/dγ is upper bounded by a finite quantity. Hence x(γ) converges.
When dxkγ /dγ is positive, the discussion is similar. This confirms convergence.
Now we prove the limiting vector limγ→+∞ x(γ) is x . Let x(∞) = limγ→+∞ x(γ). Let
f (x, γ) denote the objective function in (4.9). If x = x(∞) and by Assumption 1 the
solution to problem (4.9) is unique, we have
1 2 3
f (x(γ), γ) ≤ f (x , γ) ≤ f (x , ∞) ≤ f (x(∞), ∞), (4.17)
where inequality 1 is true because x(γ) is the minimizer at γ, inequality 2 is true because
function x1 is always larger than ρ(x, γ) (see Figure 4.3), and inequality 3 is a special
4.12. PROOFS 101
Defining L(x) = ΠN
i=1 (xi + δ), we have
We can check that f (0) = 0, f (0) = 0, and f (ε) < 0 for |ε| < 1. Hence f (ε) → 1
implies ε → 0. Hence L(x(k+1) )/L(x(k) ) → 1 implies f (ε) → 1, which implies ε → 0,
(k+1)
xi +δ
which is equivalent to (k) → 1.
xi +δ
2
102 CHAPTER 4. COMBINED IMAGE REPRESENTATION
We only need to check that the stationary point, denoted by x(∗) , of the algorithm (RIAλ )
satisfies the first-order condition (FOC) of the optimization problem (LOλ ).
If x(∗) is a stationary point of (RIAλ ), then
⎛ (∗) (∗)
⎞
sign(x1 )/(| sign(x1 )| + δ)
⎜ .. ⎟
0 = 2ΦT (Φx(∗) − y) + λ ⎜
⎝ . ⎟,
⎠
(∗) (∗)
sign(xN )/(| sign(xN )| + δ)
(∗)
where xi denotes the ith component of x(∗) . It is easy to check that the above equality is
also the FOC of a local minimum of (LOλ ). 2
In fact, we can verify that
sign(x) sign(x)
= .
| sign(x)| + δ 1+δ
Chapter 5
Iterative Methods
This chapter discusses the algorithms we use to solve for the Newton direction (see (4.12))
in our sparse representation problem. We choose an iterative method because there is a fast
algorithm to implement the matrix-vector multiplication. Since our matrix is Hermitian
(moreover symmetric), we basically choose between CG and MINRES. We choose LSQR,
which is a variation of the CG, because our system is at least positive semidefinite and
LSQR is robust against rounding error caused by finite-precision arithmetic.
In Section 5.1, we start with an overview of the iterative methods. Section 5.2 explains
why we favor LSQR and how to apply it to our problem. Section 5.3 gives some details on
MINRES. Section 5.4 contains some discussion.
5.1 Overview
where y is the vectorized analyzed image, Φ is a (flat) matrix with each column a vectorized
basis function of a certain image analysis transform (e.g., a vectorized basis function for
the 2-D DCT or the 2-D wavelet transform), x is the coefficient vector, and ρ is a separable
convex function, ρ(x) = N
i=1 ρ̄(xi ), where ρ̄ is a convex 1-D function.
We apply a damped Newton method. To find the Newton direction, the following system
103
104 CHAPTER 5. ITERATIVE METHODS
where H(xk ) and g(xk ) are the Hessian (as in (4.8)) and the gradient (as in (4.7)) of f (x)
at xk as defined in the previous chapter.
We have the following observations:
[O1] If Φ is a concatenation of several transform matrices and each of them has a fast
algorithm, then the matrix-vector multiplications with Φ, ΦT and ΦT Φ have fast
algorithms.
[O2] There is a closed form for ρ̄ (xik ) and ρ̄ (xik ), i = 1, 2, 3, . . . , so there is a low complexity
O(N ) algorithm to generate the diagonal matrix made by ρ̄ (xik ) in the Hessian and
the vector made by ρ̄ (xik ) in the gradient.
[O3] Based on [O1] and [O2], there is a low-complexity algorithm to compute the gradient
vector g(xk ).
[O4] ΦT Φ and H(xk ) are dense matrices, but H(xk ) is positive definite for 0 ≤ γ < +∞.
We choose an iterative method to solve (5.1) because the main work involves matrix-vector
multiplication, for which we have fast algorithms, and some vector inner products. For
other methods like Cholesky factorization, because of [O4], in general there will be no
low-complexity algorithms.
Adopting a general notation, we consider solving a system of linear equations
Ax = b. (5.2)
When A is symmetric and positive definite, there are two important iterative methods: CG
and MINRES.
CG Conjugate gradient (CG) methods minimize the A-norm of the error, ek A = A−1 b −
xk , b − Axk 1/2 , where ek is the error vector at step k, xk is the solution estimate at
5.1. OVERVIEW 105
MINRES The minimum residual (MINRES) method minimize the Euclidean norm of the
residual, b − Axk , at step k.
In general, [78, Page 92], when the matrix is Hermitian, MINRES is preferred in theory
because of the following inequality relation [78, Page 94]:
rkM
rkC = / M 02 ≥ rk ,
M
1 − rk /rk−1
M
where rkC is the residual at step k of the CG method and rkM is the residual at step k in
MINRES. In words, at the same step, the Euclidean norm of the residual from CG is always
larger than the Euclidean norm of the residual from MINRES.
In practical numerical computing, CG and MINRES will not perform as well as we
predict under exact arithmetic. We choose a variation of CG—LSQR—which is proven to
be more robust in finite-precision computing [117, 116]. For least-squares problems where
A = B T B and b = B T c, applying LSQR to min Bx − c22 is analytically equivalent to
CG on the normal equations B T Bx = B T c, so in the following discussion about theoretical
result, we will only mention CG rather than LSQR.
A powerful technique in analyzing the convergence rate for both CG and MINRES is the
minimax polynomial of eigenvalues. Here we summarize the key results:
[C1] The A-norm of the error in the CG algorithm for a Hermitian and positive definite
matrix A, and the Euclidean norm of the residual in the MINRES algorithm for a
general Hermitian matrix, are minimized over the spaces
e0 + span{Ae0 , A2 e0 , . . . , Ak e0 }
and
r0 + span{Ar0 , A2 r0 , . . . , Ak r0 },
respectively. These two spaces are called Krylov subspaces, and the corresponding
106 CHAPTER 5. ITERATIVE METHODS
[C2] At step k, the CG error vector and the MINRES residual vector can be written as
ek = pC
k (A)e0 ,
rk = pM
k (A)r0 ,
where pC M
k and pk are two polynomials with degree no higher than k that take value
1 at the origin. Moreover,
[C3] Suppose for a Hermitian and positive semidefinite matrix A, A = U ΛU T is its eigen-
decomposition, where U is an orthogonal matrix and Λ is a diagonal matrix. Suppose
Λ = diag{λ1 , . . . , λN }. We have sharp bounds for the norms of the error and the
residual in CG and MINRES:
ek A /e0 A ≤ min max |pk (λi )|, for CG; (5.3)
pk i=1,2,... ,N
rk /r0 ≤ min max |pk (λi )|, for MINRES. (5.4)
pk i=1,2,... ,N
In (5.3) and (5.4), if the eigenvalues are tightly clustered around a single point (away
from the origin), then the right-hand sides are more likely to be minimized; hence, iterative
methods tends to converge quickly. On the other hand, if the eigenvalues are widely spread,
especially if they lie on the both sides of the origin, then the values on the right-hand sides
of both inequalities are difficult to be minimized; hence, iterative methods may converge
slowly.
Note that the above results are based on exact arithmetic. In finite-precision computa-
tion, these error bounds are in general not true, because the round-off errors that are due
to finite precision may destroy assumed properties in the methods (e.g., orthogonality). For
further discussion, we refer to Chapter 4 of [78] and the references therein.
5.1. OVERVIEW 107
5.1.4 Preconditioner
Before we move into detailed discussion, we would like to point out that the content of this
section applies to both CG and MINRES.
When the original matrix A does not have a good eigenvalue distribution, a precondi-
tioner may help. In our case, we only need to consider the problem as in (5.2) for Hermitian
matrices. We consider solving the preconditioned system
Before giving a detailed discussion about preconditioners, we need to restate our setting.
Here A is the Hessian at step k: A = H(xk ). To simplify the discussion, we assume that Φ
is a concatenation of several orthogonal matrices:
Φ = [T1 , T2 , . . . , Tm ],
where
⎛ 1+(i−1)n2
⎞
ρ̄ (xk )
1 ⎜ ⎟
Di = λ ⎜ ..
. ⎟, i = 1, 2, . . . , m,
2 ⎝ ⎠
n2 +(i−1)n2
ρ̄ (xk )
are diagonal matrices. To simplify (5.6), we consider A left and right multiplied by a block
108 CHAPTER 5. ITERATIVE METHODS
The main result is that we found the preconditioner 1 and preconditioner 2 are not “opti-
mal”. Here “optimal” means that the matrix-vector multiplication with the matrix asso-
ciated with the preconditioner L−H should still have fast algorithms, and these fast algo-
rithms should be based on fast algorithms for matrix-vector multiplication for the matrices
T1 , T2 , . . . , Tm . So the block diagonal preconditioner is the only one we are going to use.
We describe the results about the block diagonal preconditioner here, and postpone the
discussion about the preconditioners in case 1 and 2 to Section 5.4.
The optimal block diagonal preconditioner for A1 is
⎛ ⎞
(I + S1 )−1/2
⎜ ⎟
⎜ (I + S2 )−1/2 ⎟
⎜ ⎟
⎜ .. ⎟. (5.7)
⎜ . ⎟
⎝ ⎠
(I + Sm )−1/2
A striking result is due to Demmel [78, Page 168, Theorem 10.5.3]. The key idea is that
5.2. LSQR 109
among all the block diagonal preconditioners, the one that takes the Cholesky factorizer of
the diagonal submatrix as its diagonal submatrix is nearly optimal. Here “nearly” means
that the resulting condition number cannot be larger than m times the best achievable con-
dition number by using a block diagonal preconditioner. Obviously, we have fast algorithms
to multiply with matrix (I + Si )−1/2 , i = 1, 2, . . . , m.
5.2 LSQR
LSQR [117, 116] solves the following two least-squares (LS) problems, depending on whether
the damping parameter α is zero or not.
[N] When the damping parameter is equal to zero (α = 0), solve Ax = b or minimizex b −
Ax2 .
where d is the Newton direction that we want to solve for and the remaining variables are
defined in the previous chapter. The above equation is equivalent to solving an LS problem:
,! " ! ",
, Φ y − Φxk ,
, ,
(LS) : minimize , d− , ,
d , D(xk ) −D (xk )g(xk ) ,
−1
2
110 CHAPTER 5. ITERATIVE METHODS
where
⎛ ⎞1/2
ρ̄ (xk1 )
. ⎜ ⎟
D(xk ) = λ/2 ⎜
⎝
..
. ⎟
⎠ ,
ρ̄ (xkN )
and
⎛ ⎞
ρ̄ (xk1 )
λ⎜ .. ⎟
g(xk ) = ⎜ ⎟.
2⎝ . ⎠
ρ̄ (xN ) k
Note that this is the [N] case. Recall ρ̄ and ρ̄ are defined in (4.4) and (4.5).
We can transfer problem (LS) into a damped LS problem by solving for a shifted variable
¯ δ d¯ = D(xk )d + D−1 (xk )g(xk ). The problem (LS) becomes
d:
,! " ! ",
, ΦD−1 (x )δ − − −2 (x )g(x )) ,
, k y Φ(xk D k k ,
(dLS) : minimize , d¯ − , .
d¯ , δI 0 ,
2
This is the [R] case. This method is also called diagonal preconditioning. Potentially, it
will turn the original LS problem into one with clustered singular values, so that LSQR
may take fewer iterations. LSQR also works with shorter vectors (as it handles the δI term
implicitly).
A formal description of LSQR is given in [117, page 50]. We list it here for the convenience
of readers.
1. Initialize.
β1 u1 = b, α1 v1 = AT u1 , w1 = v1 , x0 = 0, φ¯1 = β1 , ρ¯1 = α1 .
2. For i = 1, 2, 3, . . .
(a) Continue the bidiagonalization.
i. βi+1 ui+1 = Avi − αi ui
ii. αi+1 vi+1 = AT ui+1 − βi+1 vi .
(b) Construct and apply next orthogonal transformation.
i. ρi = (ρ̄i 2 + βi+1
2 )1/2
5.2.4 Discussion
There are two possible dangers in the previous approaches (LS) and (dLS). They are both
caused by the existence of large entries in the vector D−1 (xk )g(xk ). The first danger
occurs in (LS), when D−1 (xk )g(xk ) is large, the right-hand side is large even though the
elements of d will be small as Newton’s method converges. Converting (5.1) to (LS) is a
112 CHAPTER 5. ITERATIVE METHODS
somewhat unstable transformation. The second danger occurs in (dLS): when an entry of
D−1 (xk )g(xk ) is large, because d = D(xk )−1 d¯− D−2 (xk )g(xk ), there must be “catastrophic
cancellation” in some of the elements of d as the iterative method converges.
As we know, the jth entry of D−1 (xk )g(xk ) is
. . 1 2 3
−γ|xkj |
k
− 1 /e−γ|xj |/2
k
k
λ/2ρ̄ (xj )/ ρ̄ (xj ) = k
λ/2 sign(xj ) e
γ
. 1 2 γ k γ k
3
− 2 |xj | |xj |
= k
λ/2 sign(xj ) e −e 2 .
γ
Since γ usually is large in our problem, when |xkj | is significantly nonzero, the jth entry
of D−1 (xk )g(xk ) is big. It is unavoidable to have significantly nonzero |xkj |. But hopefully
the images that we consider have intrinsic sparse decompositions, hence the proportion of
significantly nonzero entries in |xkj | is small and the numerical Newton’s direction is still
accurate enough, so that our iterative algorithm will still converges to the minimum.
Solving either (LS) or (dLS) gives an approach to finding a Newton’s direction. It is
hard to tell from theory which method works better. In our numerical experiments of image
decompositions, we choose to solve (LS).
5.3 MINRES
In this section, we describe the MINRES algorithm and its preconditioned version because
MINRES has a nice theoretical property in reducing the residual norm monotonically, and
it avoids the large numbers involved in transforming (5.1) to (LS) or (dLS). We did not
implement MINRES, but a comparison with LSQR on our problem might be an interesting
research topic.
A formal description of MINRES is given on page 44 of [78]. We list it for the convenience
of the readers.
Algorithm: MINRES
2. For k = 1, 2, . . . ,
If k > 1, then
T (k − 1, k) ck−1 sk−1 T (k − 1, k)
← .
T (k, k) −s̄k−1 ck−1 T (k, k)
Set
(c) Compute pk−1 = [qk −T (k −1, k)pk−2 −T (k −2, k)pk−3 ]/T (k, k), where
undefined terms are zeros for k < 2.
MINRES can be used with a block diagonal preconditioner as specified in (5.7). A precon-
ditioned MINRES is called PMINRES. Going back to (5.5), we have
⎛ ⎞⎛ ⎞
(I + S1 )−1/2 T1
⎜ ⎟⎜ ⎟
1 ⎜
⎜ (I + S2 )−1/2 ⎟⎜
⎟⎜ T2 ⎟
⎟
L−1 =√ ⎜ .. ⎟⎜ .. ⎟
2⎜
⎝ . ⎟⎜
⎠⎝ . ⎟
⎠
(I + Sm )−1/2 Tm
And then
H
M = LL
⎛ ⎞
I + D1
⎜ ⎟
⎜ I + D2 ⎟
⎜ ⎟
= 2⎜ .. ⎟.
⎜ . ⎟
⎝ ⎠
I + Dm
Algorithm: PMINRES
2. For k = 1, 2, . . . ,
If k > 1, then
T (k − 1, k) ck−1 sk−1 T (k − 1, k)
← .
T (k, k) −s̄k−1 ck−1 T (k, k)
Set
(c) Compute pk−1 = [wk −T (k−1, k)pk−2 −T (k−2, k)pk−3 ]/T (k, k), where
undefined terms are zeros for k < 2.
(d) Set xk = xk−1 + ak−1 pk−1 , where ak−1 = βξ(k).
5.4 Discussion
1 is not op-
We argue that a preconditioner based on a complete Cholesky factorization of A
timal, because the resulting preconditioner may not have fast matrix-vector multiplication.
1 = LLT , here L denotes a lower triangular
Suppose we have the Cholesky factorization A
matrix with elements (could be block matrices) aij :
⎛ ⎞
a11
⎜ ⎟
⎜ a21 a22 ⎟
⎜ ⎟
L=⎜ ⎟.
⎜ a31 a32 a33 ⎟
⎝ ⎠
..
··· .
a11 aT11 = I + S1 ,
a11 aT12 = I,
···
a22 aT22 + a21 aT21 = I + S2 ,
···
Obviously,
a11 = (I + S1 )1/2 ,
a12 = (I + S1 )−1/2 ,
and a22 = (I + S2 − (I + S1 )−1 )−1/2 .
a22 is from the Cholesky factorization of the Schur complement of the first 2 × 2 block.
5.4. DISCUSSION 117
From
and there are fast algorithms to implement matrix-vector multiplication with matrix T1 ,
T1T , (I + D1 )−1/2 and (I + D1 )1/2 , so there are fast algorithms to multiply with a11 and a12 .
But for a22 , none of T1 , inverse of T1 , T2 and inverse of T2 can simultaneously diagonalize
I + S2 and (I + S1 )−1 . Hence there is no trivial fast algorithm to multiply with matrix a22 .
In general, a complete Cholesky factorization will destroy the structure of matrices from
which we can have fast algorithms. The fast algorithms of matrix-vector multiplication are
so vital for solving large-scale problems with iterative methods (and an intrinsic property of
our problem is that the size of data is huge) that we do not want to sacrifice the existence
of fast algorithms. Preconditioning may reduce the number of iterations, but the amount
of computation within each iteration is increased significantly, so overall, the total amount
of computing may increase. Because of this philosophy, we stop plodding in the direction
of complete Cholesky factorization.
The idea of a sparse approximate inverse (SAI) is that if we can find an approximate
eigendecomposition of the inverse matrix, then we can use it to precondition the linear
system. More precisely, suppose Z is an orthonormal matrix and at the same time the
1
columns of Z, denoted by zi , i = 1, 2, . . . , N , are A-conjugate orthogonal to each other:
Z = [z1 , z2 , . . . , zN ], and 1 T = D,
Z AZ
1 = Z T DZ and A
where D is a diagonal matrix. We have A 1−1 = Z T D−1 Z. If we can find
1 −1/2 Z)T ≈ I, so (D−1/2 Z) is a good preconditioner.
such a Z, then (D−1/2 Z)A(D
We now consider its block matrix analogue. Actually in the previous subsection, we
have already argued that when we have the block matrix version of a preconditioner, each
element of the block matrix should be able to be associated with a fast algorithm that is
based on matrix multiplication with matrices T1 , T2 , . . . , Tm and matrix-vector multiplica-
tion with diagonal matrices. Unfortunately, a preconditioner derived by using the idea of
118 CHAPTER 5. ITERATIVE METHODS
can be factorized as
1= λ11 λ12 D1 λT11 λT21
A ,
λ21 λ22 D2 λT12 λT22
where matrix
λ11 λ12
λ21 λ22
Hence
(S2 + I − (S1 + I)−1 )−1 = λ21 D1−1 λT21 + λ22 D2−1 λT22 .
So if there were fast algorithms to multiply with matrices λ21 , λ22 and their transposes,
then there would be a fast algorithm to multiply with matrix (S2 + I − (S1 + I)−1 )−1 .
But we know that in general, there is no fast algorithm to multiply with this matrix (see
also the previous subsection). So in general, such a eigendecomposition will not give a
block preconditioner that has fast algorithms to multiply with its block elements. Hence
we proved that at least in the 2 × 2 case, an eigendecomposition of inverse approach is not
favored.
Chapter 6
Simulations
Section 6.1 describes the dictionary that we use. Section 6.2 describes our testing images.
Section 6.3 discusses the decompositions based on our approach and its implication. Sec-
tion 6.4 discusses the decay of amplitudes of coefficients and how it reflects the sparsity
in representation. Section 6.5 reports a comparison with Matching Pursuit. Section 6.6
summarizes the computing time. Section 6.7 describes the forthcoming software package
that is used for this project. Finally, Section 6.8 talks about some related efforts.
6.1 Dictionary
The dictionary we choose is a combination of an orthonormal 2-D wavelet basis and a set
of edgelet-like features.
2-D wavelets are tensor products of two 1-D wavelets. We choose a type of 1-D wavelets
that have a minimum size support for a given number of vanishing moments but are as
symmetrical as possible. This class of 1-D wavelets is called “Symmlets” in WaveLab [42].
We choose the Symmlets with 8 vanishing moments and size of the support being 16. An
illustration of some of these 2-D wavelets is in Figure 3.6.
Our “edgelet dictionary” is in fact a collection of edgelet features. See the discussion
of Sections 3.3.1–3.3.3. In Appendix B, we define a collection of linear functionals λ̃e [x]
operating on x belonging to the space of N × N images. These linear functionals are
associated with the evaluation of an approximate Radon transform as described in Appendix
B. In effect, the Riesz representers of these linear functionals, {ψ̃e (k1 , k2 ) : 0 ≤ k1 , k2 < N },
119
120 CHAPTER 6. SIMULATIONS
gives a collection of “thickened edgelets”, i.e., line segments digitally sampled and of thick-
ness a few pixels wide. Our “edgelet dictionary” is precisely this collection of representers.
We will call this an edgelet dictionary even though edgelets have been previously defined
in Section 3.3.1 as line segments in R2 rather than vectors {ψ̃e (k1 , k2 ) : 0 ≤ k1 , k2 < N };
we hope this abuse of terminology will not confuse most readers. An illustration of some of
these representers can be found in Figure B.3. An obvious drawback of this set of features
is that they have roughly the same width. This property stops the set from being tight and
makes it incapable of representing fine scale image components. Developing a similar set of
linear features with various width will be an interesting research topic. See Section 7.2.
The two sets are chosen because they possess some interesting phenomena in the images
and there are fast algorithms to carry out their discrete version transforms. At this stage,
we are mainly interested in point singularities and linear singularities in an image.
6.2 Images
• Pentagon: a pentagon;
[1] They have the features that we want to test with our dictionary. Our dictionary is
made by 2-D wavelets and edgelet-like features. 2-D wavelets resemble point singu-
larities and edgelet-like features resemble linear singularities. These features are the
key image components that we are interested in, and our approach should find them.
6.3. DECOMPOSITION 121
10 5 5 5
20 10 10 10
30 15 15 15
40 20 20 20
50 25 25 25
60 30 30 30
10 20 30 40 50 60 5 10 15 20 25 30 5 10 15 20 25 30 5 10 15 20 25 30
Car: Lichtenstein Pentagon Overlapped singularities Separate singularities
[2] They are simple but sufficient to examine whether our basic assumption is true—
that different transforms will automatically represent the corresponding features in a
sparse image decomposition.
6.3 Decomposition
A way to test whether this approach works is to see how it decomposes the image into
parts associated with included transforms. Our approach presumably provides a global
sparse representation. Because of the global sparsity, if we reconstruct part of an image by
using only coefficients associated with a certain transform, then this partial reconstruction
should be a superposition of a few atoms (from the dictionary made by representers of the
transform). So we expect to see features attributable to the associated transforms in these
partial reconstructions.
We consider decomposing an image into two parts—wavelet part and edgelet part. In
principle, we expect to see points and patches in the wavelet part and lines in the edgelet
part. Figure 6.2 shows decompositions of the four testing images. The first row is for Car.
The second, third and fourth row are for Pentagon, Overlap and Separate respectively.
In each row, from the left, the first squared sub-image is the original, the second is the
wavelet part of the image, the third is the edgelet part of the image, and the last is the
superposition of the wavelet part and the edgelet part. With appropriate parameter λ, the
122 CHAPTER 6. SIMULATIONS
1. Overall, we observe the wavelet parts possess features resembling points (fine scale
wavelets) and patches (coarse scale scaling functions), while the edgelet parts possess
features resembling lines. This is most obvious in Car and least obvious in Pentagon.
The reason could be that Pentagon does not contain many linear features. (Bound-
aries are not lines.)
2. We observe some artifacts in the decompositions. For example, in the wavelet part of
Car, we see a lot of fine scale features that are not similar to points, but are similar to
line segments. This implies that they are made by many small wavelets. The reason
for this is the intrinsic disadvantage of the edgelet-like transform that we have used.
As in Figure B.3, the representers of our edgelet-like transform have a fixed width.
This prevents it from efficiently representing narrow features. The same phenomena
emerge in Overlap and Separate. A way to overcome it is to develop a transform
whose representers have not only various locations, lengths and orientations, but also
various widths. This will be an interesting topic of future research.
1. The 2-D DCT, especially the one localized to 8 × 8 blocks, has been applied in an
important industry standard for image compression and transmission known as JPEG.
2. The 2-D wavelet transform is a modern alternative to 2-D DCT. In some developing
industry standards (e.g., draft standards for JPEG-2000), the 2-D wavelet transform
has been adopted as an option. It has been proven to be more efficient than 2-D DCT
in many cases.
6.4. DECAY OF COEFFICIENTS 123
10
20
30
40
50
60
50 100 150 200 250
Original; Wavelet+Edgelet=both
5
10
15
20
25
30
20 40 60 80 100 120
Original; Wavelet+Edgelet=both
5
10
15
20
25
30
20 40 60 80 100 120
Original; Wavelet+Edgelet=both
5
10
15
20
25
30
20 40 60 80 100 120
Original; Wavelet+Edgelet=both
0 0
10 10
−1
10
−1
log (|c| ) 10
log10(|c|(i))
10 (i)
−2
10
−2
10
−3
10
−3 −4
10 10
0 100 200 300 400 500 0 100 200 300 400 500
i i
0 0
10 10
−1 −1
10 10
log (|c| )
log10(|c|(i))
10 (i)
−2 −2
10 10
−3 −3
10 10
−4 −4
10 10
0 100 200 300 400 500 0 100 200 300 400 500
i i
The goal of this thesis is to find a sparse representation of an image. We aim to see if our
combined approach will lead to a sparser representation than these existing approaches.
Figure 6.3 shows the decay of the amplitudes of coefficients from three approaches on
four images. On each plot, the horizontal axis is the order index of the amplitudes of
the coefficients and the vertical axis is the logarithm of the amplitude (in base 10). The
amplitudes are sorted from largest to smallest. The dashed lines “- - -” illustrate the
DCT-only approach; the solid lines denote our combined (wavelets+edgelet-like features)
approach; the dash and dotted lines “- · - · -” illustrate the wavelet-only approach. From
left to right and upper to lower, these plots give results for Car, Pentagon, Overlap and
Separate.
We have the following observations:
1. Our combined approach tends to give the sparsest representation (in the sense of the
fastest decay of the amplitudes) among all the three approaches. This is particularly
clear in the case of Car and Separate. But in some cases the wavelet-only approach
provides very competitive results—for example, for Pentagon and Overlap. Actually,
in Pentagon and Overlap, we observe that the curve associated with the wavelet-only
6.5. COMPARISON WITH MATCHING PURSUIT 125
approach ultimately falls well below the curve associated with the combined approach.
(This seems to imply that the wavelet-only approach gives a sparser asymptotic rep-
resentation.)
2. The DCT only approach always gives the least sparse coefficients. As we have men-
tioned, the DCT is good for images with homogeneous components. Unfortunately,
in our examples, none of them seems to have a high proportion of homogeneous com-
ponents. This may explain why here DCT is far from optimal.
3. For more “realistic” images, as in the case of Car, the wavelet-only approach works
only as well as the DCT-only approach, but our combined approach is significantly
better. This may imply that our approach is better-suited for natural images than
existing methods. Of course more careful study and extensive experiments are required
to verify this statement.
Our approach is a global approach, in the sense that we minimize a global objective func-
tion. Our method is computationally expensive. A potentially cheaper method is Matching
Pursuit (MP) [102]. MP is a greedy algorithm. In a Hilbert space, MP at every step picks
up the atom that is the most correlated with the residual at that step. But MP runs the
danger of being trapped by unfortunate choices at early steps into badly suboptimal de-
compositions [27, 40, 25]. We examine the decay of the amplitudes of coefficients for both
MP and our approach.
126 CHAPTER 6. SIMULATIONS
0 0
10 10
−1
10
−1
10
−2
10
greedy
−3
10
−2
greedy
10
global
global −4
10
−3 −5
10 10
0 200 400 600 800 1000 0 200 400 600 800 1000
0 0
10 10
−1 −1
10 10
−2 −2
10 10
10
−3 greedy −3
10 greedy
global
−4 −4
global
10 10
−5 −5
10 10
0 200 400 600 800 1000 0 200 400 600 800 1000
Figure 6.4 shows the decay of the amplitudes of the coefficients from both of the two
approaches. In these plots, the solid lines always correspond to the MP, and these dashed
lines correspond to our (minimizing the 1 norm) approach. From upper row to lower row,
left to right, the plots are for Car, Pentagon, Overlap and Separate.
We have the following observations:
1. In all four cases, our global approach always provides a sparser representation than
MP does: the decay of the amplitudes of coefficients for our approach is faster than for
MP. This validates our belief that a global optimization scheme should often provide
a sparser representation than one from a greedy algorithm.
2. If we adopt the idea that the decay of the amplitudes of coefficients at the beginning is
important, then our global approach shows the largest advantage in the example of a
realistic image (case Car). This is promising because it may imply that our approach
is more adaptable to real images.
3. We observe that our approach achieves a slightly greater advantage for Overlap than
for Separate. As we know, MP is good for separated features but not for overlapped
ones. This belief is confirmed here.
6.6. SUMMARY OF COMPUTATIONAL EXPERIMENTS 127
We ran our experiments on an SGI Power Challenger server with 196 MHz CPU. In the
algorithm in Chapter 4, we start with γ = 102 and stop with γ = 105 . Table 6.1 gives
the number of LSQR iterations and the execution time of one iteration in each case. Each
LSQR iteration includes one analysis transform and one synthesis transform of an image.
Since the number of LSQR iterations is machine independent, it is a good criterion for
comparison. The last column lists the total execution time.
In the above simulations, the tolerance parameter for LSQR is chosen to be 10−10 [117].
The minimum tolerance for Newton’s method is chosen to be 10−7 . (The Newton iteration
will stop if the norm of the gradient vector is less than 10−7 .)
6.7 Software
• MEXSource—CMEX source code for the direct edgelet transform, the fast approximate
edgelet transform, 2-D DCT, and more.
This package depends heavily on WaveLab [42] and we may integrate it as part of WaveLab.
1
An unofficial version: We started this project about two years ago. As usual, it was mixed with countless
useful and useless diversions and failures. During the past two years, I have written tens of thousands of
lines of C code and thousands of Matlab functions. There were many exciting and sleep-deprived nights.
Unfortunately, only a small proportion of the effort became the work that is presented in this thesis.
Chapter 7
Future Work
In the future, there are three promising directions that we should explore. Section 7.1 dis-
cusses more experiments that we should try but, due to the time constraints, we have not
yet performed. Section 7.2 discusses how to apply the filtering idea from multiresolution
analysis and the idea from monoscale orthonormal ridgelets to design a system with rep-
resenters having various widths. Section 7.3 talks about how to use the block coordinate
relaxation to accelerate our iterative algorithm.
7.1 Experiments
So far we have experimented on four images. We are restricted by the fact that each
experiment takes a long time to compute. To see if our approach gives sparse atomic
decomposition in more general cases, we should test it on more images. Some sets of images
that may be useful:
• Some images favored by the image processing community, for example, “barbara”,
“lenna”, etc. [2];
129
130 CHAPTER 7. FUTURE WORK
These images possess a variety of features. The results of a sparse decomposition can be
used to determine which class of features is dominant in an image, and then to tell which
transform is well-suited to processing the image.
We are going to experiment on different dictionaries. We have experimented on a dic-
tionary that is a combination of 2-D wavelets and edgelet-like features. Some other possible
combinations are:
• {2-D DCT + 2-D wavelets + edgelets}. This dictionary contains elements for ho-
mogeneous image components, point singularities and linear singularities. It should
provides sparser decomposition, but will increase the computational cost.
• {2-D DCT + 2-D wavelets + curvelets}. It is similar to the idea for the previous
dictionary, but with a better-desiged curvelets set in place of the edgelets set, it may
lead to an improvement in the computational efficiency.
• {2-D DCT + curvelets}. Comparing with the previous result, we may be able to
tell how important the wavelet components are in the image, by knowing how many
wavelets we need in sparsely representing the desired image. We can explore the same
idea for 2-D DCT and curvelets, respectively.
There are some open questions. We hope to gain more insights (and hopefully answer
the questions) via more computational experiments and theoretical analysis. Some of these
open questions are:
• We want to explore the limit of our approach in finding a sparse atomic decomposition.
For example, if the desired image can be made by a few atoms from a dictionary, will
our approach find the sparsest atomic decomposition in the dictionary?
• For most of natural images (e.g., those images we have seen in our ordinary life), can
they be represented by a few components from a dictionary made by 2-D DCT, 2-D
wavelets, edgelets and curvelets? If not, what are other transforms we should bring
in (or develop)?
• We think our sparse image representation approach can be used as a tool to preprocess
a class of images. Based on the sparse representation results, we can decide which
image transform should be chosen for the specific class of images. We are going to
explore this idea. This gives a way of doing method selection and it has applications
in image search, biological image analysis, etc.
7.2. MODIFYING EDGELET DICTIONARY 131
Suppose the matrix Φ can be split into several orthogonal and square submatrices, Φ =
[Φ1 , Φ2 , . . . , Φk ]; and x can be split into some subvectors accordingly, x = (xT1 , xT2 , . . . , xTk )T .
Problem (B1) is equivalent to
k
k
(B2) minimize y − Φi xi 22 + λ xi 1 .
x1 ,x2 ,... ,xk
i=1 i=1
1i of (B2) for xi is
When x1 , . . . , xi−1 , xi+1 , . . . , xk are fixed but not xi , the minimizer x
⎛ ⎛ ⎞⎞
⎜ ⎜
k
⎟⎟
1i = ηλ/2 ⎝ΦTi ⎝y −
x Φl xl ⎠⎠ ,
l=1
l=i
We may iteratively solve problem (B2) and at each iteration, we rotate the soft-thresholding
scheme through all subsystems. (Each subsystem is associated with a pair (Φi , xi ).) This
method is called block coordinate relaxation (BCR) in [124, 125]. Bruce, Sardy and Tseng
[124] report that in their experiments, BCR is faster than the interior point method proposed
by Chen, Donoho and Saunders [27]. They also note that if some subsystems are not
orthogonal, then BCR does not apply.
Motivated by BCR, we may split our original optimization problem into several sub-
problems. (Recall that our matrix Φ is also a combination of submatrices associated with
some image transforms.) We can develop another iterative method. In each iteration, we
solve each subproblem one by one by assuming that the coefficients associated with other
subsystems are fixed. If a subproblem corresponds to an orthogonal matrix, then the solu-
tion is simply a result of soft-thresholding of the analysis transform of the residual image. If
a subproblem does not correspond to an orthogonal matrix, moreover, if it corresponds to a
7.3. ACCELERATING THE ITERATIVE ALGORITHM 133
overcomplete system, then we use our approach (that uses Newton method and LSQR) to
solve it. Compared with the original problem, each subproblem has a smaller size, so hope-
fully this splitting approach will give us faster convergence than our previously presented
global approach (that minimizes the objective function as whole).
Some experiments with BCR and its comparison with our approach is an interesting
topic for future research.
134 CHAPTER 7. FUTURE WORK
Appendix A
This chapter is about the implementation of edgelet transform described in [50]. We give
a review of edgelets in Section A.1, then some examples in Section A.2, and finally some
details in Section A.3.
A.1 Introduction
The edgelet [50] transform was developed to represent needle-like features in images. Edgelets
are 2-D objects taking various scales, locations and orientations. If we consider an image
as a function on a unit square [0, 1] × [0, 1], an edgelet system is constructed as follows:
[E1] Partition the unit square into dyadic sub-squares, so that the sides of these squares
take values at 1/2, 1/4, 1/8, . . . .
[E2] On each dyadic subsquare, put equally-spaced vertices on the boundary, starting from
corners. We require each side equally partitioned by these vertices, and we generally
assume that there are dyadic and integral number of vertices on each side, so the
distance between two neighbor vertices should be a dyadic value too.
[E3] For a line segment that connects two vertices as in [E2], if it does not coincide with a
boundary, then it is called an edgelet.
135
136 APPENDIX A. DIRECT EDGELET TRANSFORM
A.2 Examples
Before we present details, let’s first look at some examples. The key idea of developing this
transform is hoping that if the original image is made by a few needle-like components,
then this transform will give a small number of significant coefficients, and the rest of the
coefficients will be relatively small. Moreover, if we apply the adjoint transform to the co-
efficients selected by keeping only these with significant amplitudes, then the reconstructed
image should be close to the original.
To test the above idea, we select four images:
[Huo] was selected because it is made by a few lines, so it is an ideal testing image.
[Sticky] has a patch in the head. We apply an edge filter to this image before we apply the
A.2. EXAMPLES 137
edgelet transform. [WoodGrain] is a natural image, but it has significant linear features in
it. [Lenna] is a standard testing image in image processing. As for [Sticky], we apply an
edge filter to [Lenna] before doing the edgelet transform. The images in the above list are
roughly ranked by the abundance (of course this could be subjective) of linear features.
Figure A.1, A.2, A.3 and A.4 show the numerical results. From the reconstructions
based on partial coefficients—Figure A.1 (b), (c), (e) and (f), Figure A.2 (d), (e) and (f),
Figure A.3 (b), (c), (e) and (f), Figure A.4 (d), (e) and (f)—we see that the significant
edgelet coefficients capture the linear features of the images. The following table shows the
percentages of the coefficients being used in the reconstructions. Note all the percentages
in the above table are small, say, less than 6%. The larger the percentage is, the better the
reconstruction captures the linear features in the images.
Here we design an edge filter. Let A represent the original image. Define 2-D filters D1 , D2
and D3 as
− + − − − +
D1 = , D2 = , D3 = .
− + + + + −
Let “” denote 2-D convolution. The notation [·]·2 means square each element of the matrix
in the brackets. Let “+” be an elementwise addition operator. The edge filtered image of
A is defined as
As we have explained, for the [Sticky] and [Lenna] image, we first apply an edge filter.
138 APPENDIX A. DIRECT EDGELET TRANSFORM
10 10 10
20 20 20
30 30 30
40 40 40
50 50 50
60 60 60
20 40 60 20 40 60 20 40 60
0 20 40 60 20 40 60
0 5000 10000
Figure A.1: Edgelet transform of the Chinese character “Huo”: (a) is the original; (d) is
the sorted coefficients; (b), (c), (e) and (f) are reconstructions based on the largest 200,
400, 800, 1600 coefficients, respectively.
The reason to do this is that both of them show some patchy patterns in the original images.
The filtered images show more obvious edge features. See Figure A.2 (b) and Figure A.4
(b).
A.3 Details
.
A direct way of computing edgelet coefficients is explained. It is direct in the sense that
every single edgelet coefficient is computed by a direct method; the coefficient is a weighted
sum of the pixel values (intensities) that the edgelet trespasses. For an N × N image, the
order of the complexity of the direct method is O(N 3 log N ).
Section A.3.1 gives the definition of the direct edgelet transform. Section A.3.2 calculates
the cardinality of the edgelet system. Section A.3.3 specifies the ordering. Section A.3.4
A.3. DETAILS 139
20 20 4
40 40
3
60 60
80 80 2
100 100
1
120 120
20 40 60 80 100 120 20 40 60 80 100 120 0
0 5000 10000
(d) Edgelets assoc. with largest−100 coeff. (e) Largest−300 coeff. (f) Largest−500 coeff.
20 20 20
40 40 40
60 60 60
80 80 80
100 100 100
120 120 120
20 40 60 80 100 120 20 40 60 80 100 120 20 40 60 80 100 120
Figure A.2: Edgelet transform of the sticky image: (a) is the original; (b) is the filtered
image; (c) is the sorted coefficients; (d), (e) and (f) are the reconstructions based on the
largest 100, 300 and 500 coefficients, respectively.
explains how to compute a single edgelet coefficient. Section A.3.5 shows how to do the
adjoint transform.
A.3.1 Definition
(a) Wood Grain Image (b) Largest−10000 coeff. (log) (c) Largest−20000 coeff.
210
100 100
400 400
180
500 500
170 100 200 300 400 500 100 200 300 400 500
5 10 15
5
x 10
Figure A.3: Edgelet transform of the wood grain image: (a) is the original; (d) is the sorted
coefficients; (b), (c), (e) and (f) are reconstructions based on the largest 1 × 104 , 2 × 104 ,
4 × 104 , 8 × 104 coefficients, respectively.
so when the image is a gray-scale image, I(i, j) is the gray scale at point (i, j). A 2-D
function f is defined in the following way: for i, j = 1, 2, . . . , N , if x and y fall in the
cell ((i − 1)/N, i/N ] × ((j − 1)/N, j/N ], then f (x, y) = I(i, j). Suppose e is an edgelet as
described in Section A.1. The edgelet transform of I corresponding to e is the integration
of the function f divided by the length of e:
%
ef
E(I, e) = .
length(e)
Note this integration is a weighted sum of I(i, j) where the edgelet e intersects the square
((i − 1)/N, i/N ] × ((j − 1)/N, j/N ]. The weight of I(i, j) depends on the fraction of e in
the cell ((i − 1)/N, i/N ] × ((j − 1)/N, j/N ].
A.3. DETAILS 141
100 100
200 200
300 300
4
400 400 10
500 500
100 200 300 400 500 100 200 300 400 500
1000 3000 5000
(d) Largest−1000 coeff.(log scale) (e) Largest−2000 coeff. (f) Largest−4000 coeff.
Figure A.4: Edgelet transform of Lenna image: (a) is the original Lenna image; (b) is the
filtered version; (c) is the sorted largest 5, 000 coefficients out of 428032. (d), (e) and (f)
are the reconstructions based on the largest 1000, 2000 and 4000 coefficients, respectively.
A.3.2 Cardinality
In each 2j × 2j dyadic square, every side has 1 + 2j−l vertices; hence the total number of
edgelets is
14 5
4 · (2 · 2j−l − 1) + 4 · (2j−l − 1)(3 · 2j−l − 1) = 6 · 2j−l · 2j−l − 4 · 2j−l .
2
142 APPENDIX A. DIRECT EDGELET TRANSFORM
su
#{edgelets} = 6(su − sl + 1) · 22n−2l − 4 22n−j−l
j=sl
= 6(su − sl + 1) · 2 2n−2l
− 4 · 22n−l (2−(sl −1) − 2−su ).
n
#{edgelets} = 6(n − l) · 22n−2l − 4 22n−j−l
j=l+1
= 6(n − l) · 2 2n−2l
− 4 · 22n−l (2−l − 2−n )
= 6(n − l) · 22n−2l − 22n−2l+2 + 2n−l+2
2 32 2 32 2 3
= 6(n − l) N/2l − 4 N/2l + 4 N/2l .
A.3.3 Ordering
In this subsection, we discuss how to order the edgelet coefficients. The ordering has three
layers:
There is a natural way to order scales. The scale could go from the lowest sl to the highest
su , sl ≤ su .
In the next two subsubsections, we describe how we order dyadic squares within a scale
and how we order edgelets within a dyadic square.
A.3. DETAILS 143
For scale j, we have 2n−j × 2n−j number of dyadic squares. The following table gives a
natural way to order them.
1 2 3 ··· 2n−j
2n−j + 1 2n−j + 2 2n−j + 3 ··· 2 · 2n−j
2 · 2n−j + 1 2 · 2n−j + 2 2 · 2n−j + 3 ··· 3 · 2n−j
.. .. .. .. ..
. . . . .
/ n−j 0 n−j / n−j 0 n−j / n−j 0 n−j / n−j 02
2 −1 2 +1 2 −1 2 +2 2 −1 2 + 3 ··· 2
We start with ordering the vertices on the boundary. We will see that it gives an ordering
for edgelets within this square.
For every dyadic square, the vertices on the boundary can be ordered by: starting from
the left upper corner, labeling the first vertex to the right by 1, clockwisely labeling the rest
of vertices by integers from 2 to 4 · 2j−l , where j is the scale of the dyadic square. Obviously
4 · 2j−l is the total number of vertices on this dyadic square. When n = 3 and l = 1, Figure
A.5 illustrates the ordering of vertices for a dyadic square at scale j = 3.
The ordering of vertices gives a natural ordering of edgelets. In a 2j × 2j square, there
are 4 · 2j−l vertices. Since each edgelet is determined by two vertices, let a pair (k, l) denote
the edgelet connecting vertex k and vertex l. We order the edgelets by the following two
rules:
2. Fixing first index k, let the second index l increase. When l hits the upper limit,
increase index k by one, and reset index l to the lowest possible value restricted by
the previous rule. Repeat this until none of the indices can be increased.
Combining the ordering in all three steps, we get an ordering for the edgelets. The following
system enumerates all edgelets.
144 APPENDIX A. DIRECT EDGELET TRANSFORM
16 1 2 3 4
I(1,1) I(1,2) I(1,3) I(1,4) I(1,5) I(1,6) I(1,7) I(1,8)
15 5
I(3,1) I(3,2) I(3,3) I(3,4) I(3,5) I(3,6) I(3,7) I(3,8)
14 6
I(5,1) I(5,2) I(5,3) I(5,4) I(5,5) I(5,6) I(5,7) I(5,8)
13 7
I(7,1) I(7,2) I(7,3) I(7,4) I(7,5) I(7,6) I(7,7) I(7,8)
8
12 11 10 9
Figure A.5: Vertices at scale j, for a 8 × 8 image with l = 1. The arrows shows the trend
of ordering. Integers outside are the labels of vertices.
For scale j = sl to su ,
For dyadic squares index d = 1 to 22n−2j ,
For edgelet index e = 1 to 6 · 22j−2l − 4 · 2j−l (recall l is the coarsest scale)
edgelet associated with (j, d, e);
End;
End;
End.
This section describes a direct way of computing edgelet coefficients. Basically, we compute
them one by one, taking no more than O(N ) operations each. There are O(N 2 log N ) edgelet
coefficients, so the overall complexity of a direct algorithm is no higher than O(N 3 log N ).
A.3. DETAILS 145
We have the following algorithm: for matrix I, and an edgelet with end points (x1 , y1 )
and (x2 , y2 ), where x1 , y1 , x2 and y2 are integers, we have
x
1 +1
y2
1
T(I, {x1 , y1 , x2 , y2 }) = I(i, j).
2|y2 − y1 |
i=x1 j=y1 +1
1
x2 y
1 +1
3. Case three: The edge is neither horizontal nor vertical. We must have
146 APPENDIX A. DIRECT EDGELET TRANSFORM
x1 < x2 . Let
y2 − y1
∆= .
x2 − x1
T(I, {x1 , y1 , x2 , y2 }) =
1
x2
y1 +(i−x1 )∆
ω [j, y1 + (i − x1 − 1)∆, y1 + (i − x1 )∆] I(i, j);
(y2 − y1 )
i=x1 +1 j= y1 +(i−x1 −1)∆
T(I, {x1 , y1 , x2 , y2 }) =
y1 +(i−x1 −1)∆
1
x2
ω [j, y1 + (i − x1 )∆, y1 + (i − x1 − 1)∆] I(i, j);
|y2 − y1 |
i=x1 +1 j= y1 +(i−x1 )∆
• when x = y, where x and y are the smallest integers bigger
than or equal to x and y,
(y − x) j = x,
ω [j, x, y] =
0 otherwise;
Figure A.6 gives an illustration of weights in edgelet transform. For a single edgelet
from (x1 , y1 ) to (x2 , y2 ), the weight of pixel (k, l) is equal to the length of the thick line at
the bottom of the square associated with this pixel, divided by the quantity |y2 − y1 |.
A.3. DETAILS 147
x1, y1
x ,y
2 2
In this section, we derive the adjoint edgelet transform. We start with reviewing the edgelet
transform in a symbolic way. If α is an index of the image—α = (i, j) and for an image
I, I(α) is the intensity at pixel α—and e is an index of the edgelet, the coefficient of the
edgelet transform at e is actually a weighted sum:
T(I, e) = ω(α, e)I(α).
α
Let T denote the adjoint edgelet transform. By definition, for any x ∈ RN ×N in the
image domain and any y in the edgelet transform domain, we must have
Note both sides of equation (A.1) are linear functions of y. Let y = δe , which means that
148 APPENDIX A. DIRECT EDGELET TRANSFORM
y is equal to one if and only if it is at the position corresponding to e, and zero elsewhere.
Note δe is a generalized version of the Dirac function. Equation (A.1) becomes
A.3.6 Discussion
The following observations make a fast algorithm for the edgelet transform possible:
• We can utilize inter-scale relationships. Some edgelets at coarse scales are just a linear
combination of other edgelets at a finer scale. So its coefficient is a linear combination
of other edgelet coefficients. An analogue of it is the 2-scale relationship in orthogonal
wavelets.
• We may use the special pattern of the edgelet transform. This idea is similar to the
idea in [18].
In this chapter, we present a fast algorithm to approximate the edgelet transform in discrete
cases. Note the result after the current transform is not exactly the result after a direct
edgelet transform as we presented in the previous chapter. They are close, in the sense that
we can still consider the coefficients after this transform are approximate integrations along
some line segments, but the line segments are not exactly the edgelets we described in the
previous chapter.
It is clear that in order to have a fast algorithm, it is necessary to modify the original
definition of edgelets. In many cases, there is a trade off between the simplicity or efficiency
of the algorithm and the loyalty to the original definition. The same is true here. In this
chapter, we show that we can change the system of the edgelets a little, so that a fast
(O(N 2 log N )) algorithm is feasible, and the transform still captures the linear features in
an image.
The current algorithm is based on three key foundations:
In the continuous case, the idea presented here has been extensively developed in
[48, 50, 52, 51, 49]. This approach is related to unpublished work on fast approximate
Radon transforms by Averbuch and Coifman, and to published work in the field of Syn-
thetic Aperture Radar (SAR) tomography and medical tomography. These connections to
149
150 APPENDIX B. FAST EDGELET-LIKE TRANSFORM
published work came to light only in the final stages of editing this thesis. The first Matlab
version of the algorithm being presented was coded by David Donoho. The author made
several modifications to the algorithm, and also some analysis is presented at the end of
this chapter.
The rest of the chapter is organized as following: In Section B.1, we introduce the
Fourier slice Theorem and the continuous Radon transform. Note both are for continuous
functions. Section B.2 describes in detail the main algorithm—an algorithm for fast edgelet-
like transform. The tools necessary to derive the adjoint of this transform are presented
in Section B.3. Some discussion about miscellaneous properties of the fast edgelet-like
transform are in Section B.4, including storage, computational complexity, effective region,
and ill-conditioning. Finally, we present some examples in Section B.5.
The Fourier slice theorem is the key for us to utilize the fast Fourier transform to implement
the fast Radon transform. The basic idea is that for a 2-D continuous function, if we do a
2-D Fourier transform of it, then sampling along a straight line traversing the origin, the
result is the same as the result of projecting the original function onto the same straight
line then taking 1-D Fourier transform.
We will just sketch the idea of a proof. Suppose f (x, y) is a continuous function in 2-D,
where (x, y) ∈ R2 . We use f to denote the continuous interpolation of the image I, and f
to denote a general 2-D function. Let fˆ(ξ, η) denote its 2-D Fourier transform. Then we
have
√ √
fˆ(ξ, η) = f (x, y)e−2π −1xξ −2π −1yη
e dxdy. (B.1)
y x
Taking polar coordinates in both the original domain and the Fourier domain, we have
(s)
Let fˆθ (ρ) stand for the sampling of the function fˆ along the line {(ρ cos θ, ρ sin θ) : ρ ∈
B.1. TRANSFORMS FOR 2-D CONTINUOUS FUNCTIONS 151
(s)
fˆθ (ρ) = fˆ(ρ cos θ, ρ sin θ).
Let fθ (ρ ) be the projection of f onto the line having angle θ in the original domain:
(p)
fθ (ρ )dρ
(p)
= f (x, y)dxdy. (B.2)
x cos θ +y sin θ =ρ
Later we see that this is actually the continuous Radon transform. Substituting them into
equation (B.1), we have
(s)
fˆθ (ρ) = fˆ(ρ cos θ, ρ sin θ)
= f (x, y)e−2πi(xρ cos θ+yρ sin θ) dxdy
x,y
= f (x, y)dxdye−2πiρρ
ρ x cos θ+y sin θ=ρ
fθ (ρ )e−2πiρρ dρ .
(p)
=
ρ
Note the last term is exactly the 1-D Fourier transform of the projection function fθ (ρ ).
(p)
(p)
Recalling function fθ (·) in (B.2) is defined as a projection function, we have
(p)
R(f, ρ, θ) = fθ (ρ).
152 APPENDIX B. FAST EDGELET-LIKE TRANSFORM
From the Fourier slice theorem and the above equality, we have
R(f, ρ , θ)e−2πiρρ dρ .
(s)
fˆθ (ρ) =
ρ
Thus in dimension 2, the Radon transform of f is the inverse 1-D Fourier transform of
a specially sampled (along a straight line going through the origin) 2-D Fourier transform
of the original function f . Since there are fast algorithms for the Fourier transform, we can
have fast algorithms for the Radon transform.
To some extent, the edgelet transform can be viewed as the Radon transform restricted to
a small square. There is a fast way to do the Radon transform. For an N by N image,
we can find an O(N 2 log N ) algorithm by using the fast Fourier transform. The key idea
in developing a fast approximate edgelet transform is to do a discrete version of the Radon
transform. A discrete version of the Radon transform is the algorithm we presented in this
section.
The Radon transform was originally defined for continuous functions. The Fourier slice
theorem is based on continuous functions also. Our algorithm has to be based on discrete
data, say, a matrix. A natural way to transfer a matrix to a continuous function is to
view the matrix as a sampling from a continuous function. We can calculate the analytical
Radon transform of that continuous function, then sample it to get the corresponding
discrete result.
The first question arising is how to interpolate the data. The second question is, because
the Radon transform uses polar coordinates, and a digital image is generally sampled on a
grid in the Cartesian coordinate, how do we switch the coordinates and still preserve the
fast algorithms.
Section B.2.1 gives an outline of the algorithm. Section B.2.2 describes some issues in
interpolation. Section B.2.3 is about a fast algorithm to switch from Cartesian coordinates
to polar coordinates. Section B.2.4 presents the algorithm.
B.2. DISCRETE ALGORITHM 153
B.2.1 Synopsis
We think of the Radon transform for an image as the result of the following five steps:
Outline
image
⇓ (1) Interpolate the discrete data to a continuous 2-D function.
f (x, y)
⇓ (2) Do 2-D continuous time Fourier transform.
fˆ(x, y)
⇓ (3) Switch from Cartesian to polar coordinates.
fˆ(ρ, θ)
⇓ (4) Sample at fractional frequencies according
to a grid in the polar coordinate system.
fˆ(ρi , θj )
⇓ (5) Do 1-D inverse discrete Fourier transform.
Radon transform
In subsubsection B.2.2, we describe the sampling idea associated with steps (1), (2), (4)
and (5). In subsubsection B.2.3, we describe a fast way to sample in Frequency domain.
where ρ is the interpolating kernel function. We assume ρ is equal to one at the origin and
zero at all the other integers:
notation stands for the convolution, and function δ(·) is the Dirac function at point 0.
The 2-D Fourier transform of f (x, y) is
⎛ ⎞
⎜ √ √
−1·2πξi − −1·2πηj ⎟
fˆ(ξ, η) = ⎝ I(i, j)e− e ⎠ ρ̂(ξ)ρ̂(η), (B.3)
i=1,2,... ,N ;
j=1,2,... ,N.
where ρ̂(·) is the Fourier transform of function ρ(·). Note there are two parts in fˆ(ξ, η), the
first part denoted by F (ξ, η),
√ √
F (ξ, η) = I(i, j)e− −1·2πξi − −1·2πηj
e ,
i=1,2,... ,N ;
j=1,2,... ,N.
is actually the 2-D Fourier transform of I. Note F (ξ, η) is a periodic function with period
one for both ξ and η: F (ξ + 1, η) = F (ξ, η) and F (ξ, η + 1) = F (ξ, η). If we sample ξ and η
1 2
at points N, N,... , 1, then we have the discrete Fourier transform (DFT). We know there
is an O(N log N ) algorithm to implement.
Function ρ̂(·) typically has finite support. In this paper, we choose the support to have
length equal to one, so that the support of ρ̂(ξ)ρ̂(η) forms a unit square. We did not choose
a support wider than one for some reason we will mention later.
From equation (B.3), fˆ(ξ, η) is the periodic function F (ξ, η) truncated by ρ̂(ξ)ρ̂(η). The
shape of ρ̂(ξ) and ρ̂(η) determines the property of function fˆ(ξ, η).
Section B.6.1 gives three examples of interpolation functions and some related discussion.
In this paper, we will choose the raised cosine function as our windowing function.
In this section, we describe how to transfer from Cartesian coordinates to polar coordinates.
As in the synopsis, we do the coordinate switch in the Fourier domain.
From equation (B.3), function fˆ(ξ, η) is just a multiplication of function F (ξ, η) with a
windowing function. If we know the function F (ξ, η), the function fˆ(ξ, η) is almost a direct
B.2. DISCRETE ALGORITHM 155
extension. In the following, we treat fˆ(ξ, η) as a box function (e.g., the indicator function
of the unit square).
Note the range of the index of I is changed. This is to follow the convention in DFT. By
applying the FFT, we can get the following matrix by an O(N 2 log2 N ) algorithm:
/ 0
F ( Nk , Nl ) k=0,1,... ,N −1 ,
l=0,1,... ,N −1
because for 0 ≤ k, l ≤ N − 1,
√ √
I(i + 1, j + 1)e−2π −1 N i −2π −1 Nl j
k
F ( Nk , Nl ) = e .
i=0,1,... ,N −1
j=0,1,... ,N −1
For Radon transform, instead of sampling at Cartesian grid point (k/N, l/N ), we need
to sample at polar grid points. We develop an X-interpolation approach, which is taking
samples at grid points in polar coordinate. So after X-interpolation, we get the following
N × 2N matrix:
&
/ 0
F ( Nk − 12 , s(k) + l · ∆(k) − 12 ) k=0,1,... ,N −1 ,...
l=0,1,... ,N −1
'
/ 0
F (s(k) + l · ∆(k) − 12 , NN−k − 12 ) k=0,1,... ,N −1 , (B.4)
l=0,1,... ,N −1
where
s(k) = k
N, ∆(k) = 1
N (1 − 2k
N ). (B.5)
There is a fast way to compute the matrix in (B.4), using an idea from [6].
156 APPENDIX B. FAST EDGELET-LIKE TRANSFORM
1 2 3 1 2 3 4 5 6
4 X-interpolate
5 =⇒
6
For 0 ≤ k, l ≤ N − 1, we have
k 1 1
F − , s(k) + l · ∆(k) −
N 2 2
√ √
−2π −1( N
k
− 21 )i −2π −1(s(k)+l·∆(k)− 21 )j
= I(i + 1, j + 1)e e
i=0,1,... ,N −1;
j=0,1,... ,N −1.
⎛ ⎞
√ √
⎝ I(i + 1, j + 1)e−2π −1( N − 12 )i ⎠ e−2π −1(s(k)+l·∆(k)− 21 )j
k
=
j=0,1,... ,N −1 i=0,1,... ,N −1
define as g(k, j)
√
= g(k, j)e−2π −1(s(k)+l·∆(k)− 21 )j
j=0,1,... ,N −1
√ 2 √ 2 √
−2π −1∆(k) l2 −1(j·s(k)− 21 j+∆(k) j2 )
= e · g(k, j)e−2π · eπ −1∆(k)(l−j)2
.
j=0,1,... ,N −1
convolution.
The matrix (g(k, j)) k=0,1,... ,N −1 is the 1-D Fourier transform of the image matrix I(i +
j=0,1,... ,N −1
1, j+1) i=0,1,... ,N −1 by column. We can compute matrix (g(k, j)) k=0,1,... ,N −1 with O(N 2 log N )
j=0,1,... ,N −1 j=0,1,... ,N −1
work. For fixed k, computing function value F ( Nk , s(k) + l · ∆(k)) is basically a convolution.
/ 0
To compute the kth row in matrix F ( Nk , s(k) + l · ∆(k)) k=0,1,... ,N −1 , we utilize the Toeplitz
l=0,1,... ,N −1
structure. We can compute the kth row in O(N log N ) time, and hence computing the first
half of the matrix in (B.4) has O(N 2 log N ) complexity.
B.2. DISCRETE ALGORITHM 157
For the second half of the matrix in (B.4), we have the same result:
1 N −k 1
F s(k) + l · ∆(k) − , −
2 N 2
√ √ −k
−2π −1(s(k)+l·∆(k)− 21 )i −2π −1( NN − 21 )j
= I(i + 1, j + 1)e e
i=0,1,... ,N −1;
j=0,1,... ,N −1.
⎛ ⎞
√ √
⎝
j
−1( N −1−k ⎠
= I(i + 1, j + 1)e−2π −1 N
e−2π N
− 12 )j
i=0,1,... ,N −1 j=0,1,... ,N −1
define as h(i, N − 1 − k)
√
·e−2π −1(s(k)+l·∆(k)− 2 )i
1
√ 2 ( 2 2
)3
−2π −1 i·s(k)− 12 i+∆(k) l2 + i2 − 12 (l−i)2
= h(i, N − 1 − k)e
i=0,1,... ,N −1
√ 2 √ 2 √
−2π −1∆(k) l2
= e · h(i, N − 1 − k)e−2π −1(i·s(k)− 12 i+∆(k) i2 )
· eπ −1∆(k)(l−i)2
.
i=0,1,... ,N −1
convolution.
The matrix (h(i, k)) i=0,1,... ,N −1 is an assembly of the 1-D Fourier transform of rows of the
k=0,1,... ,N −1
image matrix (I(i + 1, j + 1)) i=0,1,... ,N −1 . For the same reason, the second half of the matrix
j=0,1,... ,N −1
in (B.4) can be computed with O(N 2 log N ) complexity.
The discrete Radon transform is the 1-D inverse discrete Fourier transform of columns
of the matrix in (B.4). Obviously the complexity at this step is no higher than O(N 2 log N ),
so the overall complexity of the discrete Radon transform is O(N 2 log N ).
In order to make each column of the matrix in (B.4) be the DFT of a real sequence, for
any fixed l, 0 ≤ l ≤ N − 1, and k = 1, 2, . . . , N/2, we need to have
k 1 1 N −k 1 1
F − , s(k) + l · ∆(k) − =F − , s(N − k) + l · ∆(N − k) − ,
N 2 2 N 2 2
and
1 N −k 1 1 k 1
F s(k) + l · ∆(k) − , − = F s(N − k) + l · ∆(N − k) − , − .
2 N 2 2 N 2
It is easy to verify that for s(·) and ∆(·) in (B.5), the above equations are satisfied.
158 APPENDIX B. FAST EDGELET-LIKE TRANSFORM
B.2.4 Algorithm
Here we give the algorithm to compute the fast edgelet-like transform. Note when the
original image is N by N , our algorithm generates an N × 2N matrix. Recall that I(i, j)
denotes the image value at pixel (i, j), 1 ≤ i, j ≤ N .
1. For j = 0 to N − 1 step 1,
take DFT of (j + 1)th column:
[g(0, j), g(1, j), . . . , g(N − 1, j)]
= DF T N ([I(1, j + 1), I(2, j + 1), . . . , I(N, j + 1)]);
End;
2. For k = 0 to N − 1 step 1,
(a) pointwise multiply row k + 1, [g(k,
0), g(k, 1), .. . , g(k, N − 1)], 6
√ 2
−2π −1 j·s(k)+∆(k) j2
with complex sequence e : j = 0, 1, . . . , N − 1 ;
(c) pointwise
7 multiply the sequence with complex8 sequence
√ l2
e−2π −1∆(k) 2 : l = 0, 1, . . . , N − 1 ;
9 :
we get another row vector: F ( Nk , s(k) + l · ∆(k)) : l = 0, 1, . . . , N − 1 ;
End;
3. For l = 0 to N − 1 step 1,
taper the (l + 1)th column by the raised cosine, then take inverse 1-D
DFT;
End;
1. For i = 0 to N − 1 step 1,
For j = 0 to N − 1 step 1,
√ j
modulate: g(i, j) = I(i + 1, j + 1)e−2π −1 N
;
B.3. ADJOINT OF THE FAST TRANSFORM 159
End;
End;
2. For i = 0 to N − 1 step 1,
take 1-D DFT of the (i + 1)th row:
[g(i, 0), g(i, 1), . . . , g(i, N − 1)]
= DF T N ([I(i + 1, 1), I(i + 1, 2), . . . , I(i + 1, N )]);
End;
3. flip matrix (g(i, k))0≤i,k≤N −1 by columns:
h(i, k) = g(i, N − 1 − k);
4. For k = 0 to N − 1 step 1,
(a) pointwise multiply (k + 1)th column, {h(0, k), h(1, k), . . . , h(N −
1, k)}, 7 8
√ 2 2
3
−2π −1 i·s(k)+∆(k) i2
with complex sequence e : i = 0, 1, . . . , N − 1 ;
(b) convolve it with complex sequence
4 √ 5
eπ −1∆(k)(t) : t = 0, 1, . . . , N − 1 ;
2
(c) pointwise
7 multiply the sequence with complex8 sequence
√ l2
e−2π −1∆(k) 2 : l = 0, 1, . . . , N − 1 ;
End;
5. take transpose;
6. For l = 0 to N − 1 step 1,
taper the (l + 1)th column by the raised cosine function, then take
inverse 1-D DFT;
End;
To find the adjoint of step 2 of the algorithm in computing the first half, which is also
the adjoint of step 3 in computing the second half, we develop the following linear algebra
point of view. Let x = (x1 , x2 , . . . , xN )T be a column vector. For fixed k, after taking step
two (or three for the second half), we get column vector y = (y1 , y2 , . . . , yN )T . Vectors x
and y satisfy the following equation:
⎛ √ 2
⎞
e−2π −1∆(k) 02
⎜ ⎟
⎜ .. ⎟
y = ⎜ . ⎟
⎝ √ ⎠
(N −1)2
e−2π −1∆(k) 2
⎛ √ √ √ √ ⎞
eπ −1∆(k)·02 eπ −1∆(k)·12 eπ −1∆(k)·22 ... eπ −1∆(k)·(N −1)2
⎜ √ √ √ √ ⎟
⎜ eπ −1∆(k)·12 eπ −1∆(k)·02 eπ −1∆(k)·12 eπ −1∆(k)·(N −2)2 ⎟
⎜ ... ⎟
⎜ √ √ √ ⎟
⎜ eπ −1∆(k)·22 eπ −1∆(k)·12 eπ −1∆(k)·02 ⎟
⎜ ⎟
⎜ .. .. .. ⎟
⎜ . . . ⎟
⎝ √ √ √ ⎠
eπ −1∆(k)·(N −1)2 eπ −1∆(k)·(N −2)2 ... eπ −1∆(k)·02
⎛ √ 2
⎞
e−2π −1(0·s(k)+∆(k) 02 )
⎜ ⎟
⎜ .. ⎟
⎜ . ⎟ x.
⎝ √ ⎠
(N −1)2
e−2π −1((N −1)·s(k)+∆(k) 2
)
The adjoint transform of this step corresponds to multiplying with the adjoint of the
B.4. ANALYSIS 161
From all the above, we can derive the adjoint of the transform.
B.4 Analysis
28 times the effort of doing the 2-D fast Fourier transform for an N by N matrix.
The storage needed for the previous algorithm is proportional to N 2 .
The way we sample in the 2-D Fourier domain is actually equivalent to sampling in the
region surrounded by a dashed curve in Figure B.2. A column of the output of the fast
algorithm is an equally-spaced sample of every straight line passing through the origin
within this region. The dashed curve is made by four half circles. The reason that the
region looks like this is simple: dilation in the Fourier domain is equivalent to shrinkage in
the time domain with the same factor.
Effective region
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
The fast X-interpolation transform can be divided into three steps: if we only consider the
first half of the matrix, we have
Note the matrix associated with the fractional Fourier transform may have a large
condition number. We leave further discussion as future research.
B.5. EXAMPLES 163
B.5 Examples
For any linear transform, we can regard the coefficients as inner products of the input
signal with given basic elements. If it is an isometric transform, then these basic elements
form an orthonormal basis. Note our fast edgelet-like transform is redundant. Hence the
corresponding set of basic elements does not form a basis.
In Figure B.3, we plot some basic elements of our transform. Note we actually apply
our transform for square images with different size (scale). In the first (second, third) row,
1 1
the squared images have sides 4 ( 2 , 1) of the side of the original image. We intentionally
renormalize the basic elements so that each of them should have 2 norm roughly equal to
1.
10 10 10 10
20 20 20 20
30 30 30 30
10 20 30 10 20 30 10 20 30 10 20 30
(e) MSET (2−5) (f) MSET (2−5) (g) MSET (2−5) (h) MSET (2−5)
20 20 20 20
40 40 40 40
60 60 60 60
80 80 80 80
Since this algorithm is designed to capture the linear features in an image, it will be inter-
esting to see how it works on some artificial images made by a single linear singularity.
The first row of Figure B.4 shows some needle-like images. The second row is their
multiscale fast edgelet-like transforms. To explain the images of the coefficients, we need
to explain how we arrange the output of our algorithm. Suppose we only do the transform
for the whole image (scale-0 transform). Then as in B.2.4, the input I is mapped to a
coefficient matrix with two components [E11 , E12 ]. When we divide the image into 2 × 2
block images (scale-1 transform), each subimage is mapped to a coefficient sub-matrix with
two components:
! " ! "
I11 I12 1
E11 2
E11 1
E12 2
E12
→ .
I21 I22 1
E21 2
E21 1
E22 2
E22
B.5. EXAMPLES 165
And so on. Note for a fixed scale, the output matrix has the same number of rows, but the
number of columns expands by 2. The second row in Figure B.4 is illustrations of output
by taking four scales. So the number of rows are expanded by 4. The dark points are
associated with coefficients with large amplitudes. Obviously there are few coefficients with
significantly large amplitudes and all the rest are small. This matches our original goal to
design this transformation.
100
200
300
400
500
100 200 300 400 500
(scale=2) (scale=1) (scale=0)
4000
1600 2500
Here real means that they are from some other independent sources and are not made
intentionally for our method. In these images, the linear features are embedded, and we see
that our transform still captures them.
Figure B.5 is a wood grain image. There are many needle-like objects in the image, and
they are pretty much at the the same scale—having the same length—and along the same
direction. We do our transform at three scales (0, 1, 2 corresponding to no division of the
166 APPENDIX B. FAST EDGELET-LIKE TRANSFORM
image, divided by two, and divided by four). The upper-left image is the original. The
first row gives the columnwise maximum of the absolute values of the coefficients. We see
that there are strong patterns in the columnwise maximums; in particular, it becomes large
along the direction that is the direction of most of the needle-like features in the image.
The second row are the coefficient matrices at different scales.
It takes about 40 seconds to carry out this transform on an SGI Onyx workstation.
100
200
300
400
500
100 200 300 400 500
Another example is the Barbara image. Again, the upper-left image in Figure B.6 is the
original image. The remaining images are the coefficient matrices corresponding to scale 0
through 5. The dark area corresponds to the significant coefficients. We observe that when
the scale is increased, we have more significant coefficients. See the coefficient matrix at
scale equal to 5. This implies that when we divide the image into small squares, the linear
features become more dramatic, hence it becomes easier for a monoscale fast edgelet-like
transform to capture them. Further discussion is beyond the scope of this chapter; we will
leave it as future research.
B.6. MISCELLANEOUS 167
B.6 Miscellaneous
1, − 12 < ξ ≤ 12 ;
ρ̂(ξ) =
0, otherwise.
Figure B.7 illustrates the sinc function and its Fourier transform.
1 2
0.8
1.5
0.6
1
0.4
0.5
0.2
0
0
−0.2 −0.5
−0.4 −1
−5 0 5 −2 0 2
x ξ
Figure B.7: Sinc function in (a) and its Fourier transform—blocky function in (b).
1
2 + 12 cos(2πξ), − 12 < ξ ≤ 12 ;
ρ̂(ξ) =
0, otherwise.
168 APPENDIX B. FAST EDGELET-LIKE TRANSFORM
So ρ is
√
−1ξx sin(πx)
ρ(x) = ρ̂(ξ)e2π dξ = − .
2πx(x2 − 1)
1
ρ(0) = 1, ρ(1) = ρ(−1) = , and ρ(i) = 0, i = ±2, ±3, . . . .
2
Figure B.8 shows the raised cosine function and the corresponding ρ in the time
domain.
(a) ρ (b) ρ_hat is raised cosine
2.5
1 2
0.8
1.5
0.6
1
0.4
0.5
0.2
0
0
−0.2 −0.5
−0.4 −1
−5 0 5 −2 0 2
x ξ
Figure B.8: Raised cosine function in (b) and its counterpart in time domain in (a).
3. In the first case, there is no tapering. In the second case, the tapering rate is 100%. It
is interesting to look at something in between. Suppose p (0 < p < 12 ) is the tapering
point, which means that ρ̂ is chosen to be
⎧
⎪
⎪ 1, |ξ| < p;
⎨
π(|ξ|−p)
ρ̂(ξ) = 1
+ 12 cos , p < |ξ| < 12 ;
⎪ −p
2 1
⎪
⎩
2
0, otherwise.
In this case,
1
−p
tapering rate = 2
1 = 100(1 − 2p)%.
2
B.6. MISCELLANEOUS 169
The function ρ is
√
−1ξx
ρ(x) = ρ̂(ξ)e2π dξ
&'
p 1/2
2π(ξ − p)
= cos(2πξx)dξ + cos(2πξx) 1 + cos dξ
−p p 1 − 2p
sin(πx) + sin(2πpx) 1
= .
2πx 1 − (1 − 2p)2 x2
When p is 1/2, we get the sinc function in case one. When p is 0, we get the ρ function
in case two. When p = 14 , the Figure B.9 shows how the function ρ and function ρ̂
look.
When p is 1/4 and function ρ(·) only takes integer values, unlike the two previous
cases, the function ρ does not have finite support.
1 2
0.8
1.5
0.6
1
0.4
0.5
0.2
0
0
−0.2 −0.5
−0.4 −1
−5 0 5 −2 0 2
x ξ
Figure B.9: Fifty percent tapered window function in (b) and its counterpart in time domain
in (a).
170 APPENDIX B. FAST EDGELET-LIKE TRANSFORM
Bibliography
[3] H. K. Aghajan. Subspace Techniques for Image Understanding and Computer Vision.
PhD thesis, Stanford University, March 1995.
[6] David H. Bailey and Paul N. Swarztrauber. The fractional Fourier transform and
applications. SIAM Rev., 33(3):389–404, 1991.
[7] Richard G. Baraniuk and Douglas L. Jones. Shear madness: new orthonormal
bases and frames using chirp functions. IEEE Transactions on Signal Processing,
41(12):3543–9, December 1993.
[9] Michele Benzi, Carl D. Meyer, and Miroslav T̊uma. A sparse approximate inverse
preconditioner for the conjugate gradient method. SIAM J. Sci. Comput., 17(5):1135–
1149, 1996.
171
172 BIBLIOGRAPHY
[10] T. Berger. Rate Distortion Theory. Prentice-Hall, Englewood Cliffs, NJ, 1971.
[12] Julian Besag. Spatial interaction and the statistical analysis of lattice systems (with
discussion). J. Royal Statistical Society, Series B, Methodological, 36:192–236, 1974.
[13] G. Beylkin. Discrete radon transform. IEEE Transactions on Acoustics, Speech and
Signal Processing, ASSP-35(2):162–72, Feb. 1987.
[14] G. Beylkin. On the fast Fourier transform of functions with singularities. Applied and
Computational Harmonic Analysis, 2(4):363–81, 1995.
[15] Peter Bloomfield. Fourier Analysis of Time Series: An Introduction. John Wiley &
Sons, 1976.
[17] Ronald N. Bracewell. The Fourier Transform and Its Applications. McGraw-Hill,
1986.
[18] M. L. Brady. A fast discrete approximation algorithm for the radon transform. SIAM
J. Computing, 27(1):107–19, February 1998.
[20] C. Sidney Burrus, Ramesh A. Gopinath, and Haitao Guo. Introduction to Wavelets
and Wavelet Transforms: A Primer. Prentice Hall, 1998.
[21] Emmanuel J. Candès. Ridgelets: Theory and Applications. PhD thesis, Stanford
University, 1998.
[22] Emmanuel J. Candès. Harmonic analysis of neural networks. Applied and Computa-
tional Harmonic Analysis, 6(2):197–218, 1999.
[23] J. Capon. A probabilistic model for run-length coding of pictures. IRE Transactions
on Information Theory, pages 157–63, 1959.
BIBLIOGRAPHY 173
[24] C. Victor Chen and Hao Ling. Joint time-frequency analysis for radar signal and
image processing. IEEE Signal Processing Magazine, 16(2):81–93, March 1999.
[25] Scott Shaobing Chen. Basis Pursuit. PhD thesis, Stanford University, November
1995.
[26] Scott Shaobing Chen, David Donoho, and Michael A. Saunders. About Atomizer.
Stanford University, April 1996.
[27] Scott Shaobing Chen, David L. Donoho, and Michael A. Saunders. Atomic decompo-
sition by basis pursuit. SIAM J. Scientific Computing, 20(1):33–61, 1999. electronic.
[29] Xiru Chen, Lincheng Zhao, and Yuehua Wu. On conditions of consistency of ml1 n
estimates. Statistical Sinica, 3:9–18, 1993.
[30] Charles Chui. class notes for stat323: Wavelets and beyond, with applications. Stan-
ford University, Spring 1998.
[31] P. Concus, G. H. Golub, and G. Meurant. Block preconditioning for the conjugate
gradient method. SIAM J. Sci. Statist. Comput., 6(1):220–52, 1985.
[32] J. W. Cooley and John W. Tukey. An algorithm for the machine calculation of complex
Fourier series. Math. Comp., 19:297–301, 1965.
[33] Thomas M. Cover and J. A. Thomas. Elements of Information Theory. John Wiley
& Sons, 1991.
[35] Geoffrey Davis. Adaptive Nonlinear Approximations. PhD thesis, Courant Institute
of Mathematical Sciences, September 1994.
[36] A.H. Delaney and Y. Bresler. Globally convergent edge-preserving regularized recon-
struction: an application to limited-angle tomography. IEEE Transactions on Image
Processing, 7(2):204–21, 1998.
174 BIBLIOGRAPHY
[37] E. M. Deloraine and A. H. Reeves. The 25th anniversary of pulse code modulation.
IEEE Spectrum, pages 56–64, May 1965.
[39] R. S. Dembo and T. Steihaug. Truncated Newton algorithms for large-scale uncon-
strained optimization. Math. Programming, 26(2):190–212, 1983.
[40] R.A. DeVore and V.N. Temlyakov. Some remarks on greedy algorithms. Advances in
Computational Mathematics, 5(2-3):173–87, 1996.
[41] Y. Dodge. Statistical Data Analysis: Based on the L1 −norm and Related Methods.
North-Holland, 1987.
[42] D. L Donoho and etc. WaveLab. Stanford University, Stanford, CA, .701 edition.
[Link]
[43] David Donoho and Xiaoming Huo. Uncertainty principles and ideal atomic decompo-
sition. Working paper, June 1999.
[44] David L. Donoho. Interpolating wavelet transforms. Technical report, Stanford Uni-
versity, 1992.
[45] David L. Donoho. Smooth wavelet decompositions with blocky coefficient kernels.
Technical report, Stanford University, 1993.
[46] David L. Donoho. Unconditional bases are optimal bases for data compression and for
statistical estimation. Applied and Computational Harmonic Analysis, 1(1):100–15,
1993.
[47] David L. Donoho. Unconditional bases and bit-level compression. Applied and Com-
putational Harmonic Analysis, 3(4):388–92, 1996.
[48] David L. Donoho. Fast ridgelet transforms in dimension 2. Personal copy, April 30
1997.
[49] David L. Donoho. Digital ridgelet transform via digital polar coordinate transform.
Stanford University, 1998.
BIBLIOGRAPHY 175
[50] David L. Donoho. Fast edgelet transforms and applications. Manuscript, September
1998.
[51] David L. Donoho. Orthonormal Ridgelets and Linear Singularities. Stanford Univer-
sity, 1998.
[52] David L. Donoho. Ridge functions and orthonormal ridgelets. Stanford University,
1998. [Link]
[53] David L. Donoho. Sparse components of images and optimal atomic decom-
positions. Technical report, Stanford University, December 1998. [Link]
stat/˜donoho/Reports/1998/[Link].
[55] David L. Donoho and Stark P. B. Uncertainty principles and signal recovery. SIAM
J. Applied Mathematics, 49(3):906–31, 1989.
[56] David L. Donoho and Peter J. Huber. The notion of breakdown point. In A Festschrift
for Erich L. Lehmann, pages 157–84. Wadsworth Advanced Books and Software, 1983.
[57] David L. Donoho and B. F. Logan. Signal recovery and the large sieve. SIAM J.
Applied Mathematics, 52(2):577–91, April 1992.
[58] David L. Donoho, Stéphane Mallat, and R. von Sachs. Estimating covariance of
locally stationary processes: Consistency of best basis methods. In Proceedings of
IEEE Time-Frequency and Time-Scale Symposium, Paris, July 1996.
[59] David L. Donoho, Stéphane Mallat, and R. von Sachs. Estimating Covariance of
Locally Stationary Processes: Rates of Convergence of Best Basis Methods, February
1998.
[60] Serge Dubuc and Gilles Deslauriers, editors. Spline Functions and the Theory of
Wavelets. American Mathematical Society, 1999.
[61] P. Duhamel and C. Guillemot. Polynomial transform computation of the 2-d dct. In
ICASSP 1990, International Conference on Acoustics, Speech and Signal Processing,
volume 3, pages 1515–18, 1990.
176 BIBLIOGRAPHY
[62] P. Duhamel and H. H’Mida. New 2n DCT algorithms suitable for VLSI implemen-
tation. In Proceedings: ICASSP 87. 1987 International Conference on Acoustics,
Speech, and Signal Processing, volume 3, pages 1805–8, 1987.
[63] P. Duhamel and M. Vetterli. Improved Fourier and Hartley transform algorithms:
Application to cyclic convolution of real data. IEEE Transactions on Acoustics, Speech
and Signal Processing, ASSP-35(6):818–24, 1987.
[64] Douglas F. Elliott and K. Ramamohan Rao. Fast transforms: algorithms, analyses,
applications. Academic Press, 1982.
[66] E. Feig and S. Winograd. Fast algorithms for the discrete cosine transform. IEEE
Transactions on Signal Processing, 40(9):2174–93, September 1992.
[67] E. Feig and S. Winograd. On the multiplicative complexity of discrete cosine trans-
forms. IEEE Transactions on Information Theory, 38(4):1387–91, July 1992.
[70] A. Gersho and Robert M. Gray. Vector Quantization and Signal Compression. Kluwer
Academic Publishers, 1992.
[71] P. E. Gill, W. Murray, and M. A. Saunders. User’s Guide for SNOPT 5.3: A Fortran
Package for Large-Scale Nonlinear Programming, May 20 1998. Draft.
[72] Gene Golub, Anne Greenbaum, and Mitchell Luskin, editors. Recent Advances in
Iterative Methods, volume 60 of IMA volumes in mathematics and its applications.
Springer-Verlag, 1994. Papers from the IMA Workshop on Iterative Methods for
Sparse and Structured Problems, held in Minneapolis, Minn., Feb. 24-Mar. 1, 1992.
[73] Gene H. Golub and Dianne P. O’Leary. Some history of the conjugate gradient and
Lanczos algorithms: 1948–1976. SIAM Review, 31(1):50–102, 1989.
[74] Gene H. Golub and Charles F. Van Loan. Matrix Computations. Johns Hopkins
University Press, 3rd edition, 1996.
BIBLIOGRAPHY 177
[78] Anne Greenbaum. Iterative Methods for Solving Linear Systems, volume 17 of Fron-
tiers in Applied Mathematics. SIAM, 1997.
[81] Arthur E. Hoerl and Robert W. Kennard. Ridge regression: Biased estimation for
nonorthogonal problems. Technometrics, 12:55–67, 1970.
[83] TaiChiu Hsung, Daniel P.K. Lun, and Wan-Chi Siu. The discrete periodic radon
transform. IEEE Transactions on Signal Processing, 44(10):2651–7, October 1996.
[84] J. Huang. Quantization of Correlated Random Variables. PhD thesis, Yale University,
New Haven, CT, 1962.
[85] J.-Y. Huang and P. M. Schultheiss. Block quantization of correlated gaussian random
variables. IEEE Trans. Comm., 11:289–96, September 1963.
[86] Peter J. Huber. Fisher information and spline interpolation. Annals of Statistics,
2:1029–33, 1974.
[87] Peter J. Huber. Robust Statistical Procedures, volume 27. CBMS-NSF, 1977.
178 BIBLIOGRAPHY
[88] Peter J. Huber. Robust Statistics. Wiley Series in Pro. and Math. Sci., 1981.
[89] S. Jaggi, W.C. Karl, S. Mallat, and A.S. Willsky. High-resolution pursuit for feature
extraction. Applied and Computational Harmonic Analysis, 5(4):428–49, October
1998. [Link]
[90] R. A. Kennedy and Z. Ding. Blind adaptive equalizer for quadrature amplitude mod-
ulation communication systems based on convex cost functions. Optical Engineering,
31(6):1189–99, June 1992.
[92] B.G. Lee. A new algorithm to compute the discrete cosine transform. IEEE Trans-
actions on Acoustics, Speech and Signal Processing, ASSP-32(6):1243–5, 1984.
[93] W. Li and J.J. Swetits. The linear l1 estimator and the Huber m-estimator. SIAM J.
Optim., 8(2):457–75, May 1998.
[94] S. P. Lloyd. Least squares quantization in PCM. IEEE Trans. Inform. Th., IT-
28(2):129–37, March 1982. Originally an unpublished Bell Telephone Laboratories
tech. memo., 1957.
[95] M. Lobo and M. Fazel. Group presentation. ISL Stanford University, 1999.
[96] B.F. Logan. Properties of high-pass signals. PhD thesis, Columbia University, New
York, 1965.
[99] Stéphane Mallat. Applied mathematics meets signal processing. In Proceedings of the
International Congress of Mathematicians, Berlin, 1998.
BIBLIOGRAPHY 179
[100] Stéphane Mallat. A Wavelet Tour of Signal Processing. Academic Press, 1998.
[101] Stéphane Mallat, George Papanicolaou, and Z. Zhang. Adaptive covariance estimation
of locally stationary processes. Annals of Statistics, 26(1):1–47, February 1998.
[102] Stéphane Mallat and Zhifeng Zhang. Matching pursuits with time-frequency dictio-
naries. IEEE Transactions on Signal Processing, 41(12):3397–415, December 1993.
[103] S. Mann and S. Haykin. The chirplet transform: physical considerations. IEEE
Transactions on Signal Processing, 43(11):2745–61, November 1995.
[104] J. Max. Quantizing for minimum distortion. IRE Trans. Inform. Th., IT-6(1):7–12,
March 1960.
[105] Francois G. Meyer and Ronald R. Coifman. Brushlets: a tool for directional image
analysis and image compression. Applied and Computational Harmonic Analysis,
4:147–187, 1997.
[107] Yves Meyer. Wavelets, Algorithms & Applications. SIAM, 1993. Translated and
revised by Robert Ryan.
[108] C. Michelot and M.L. Bougeard. Duality results and proximal solutions of the Huber
M-estimator problem. Applied Mathematics & Optimization, 30:203–21, 1994.
[110] José M.F. Moura and Nikhil Balram. Recursive structure of noncausal Gauss-Markov
random fields. IEEE Transactions on Information Theory, 38(2):334–54, March 1992.
[111] José. M.F. Moura and Marcelo G.S. Bruno. DCT/DST and Gauss-Markov fields:
conditions for equivalence. IEEE Transactions on Signal Processing, 46(9):2571–4,
September 1998.
[112] Balas K. Natarajan. Sparse approximate solutions to linear systems. SIAM J. Com-
put., 24(2):227–234, 1995.
180 BIBLIOGRAPHY
[114] Henry J. Nussbaumer. Fast Fourier Transform and Convolution Algorithms. Springer-
Verlag, 1982.
[115] Kramer H. P. and Mathews M. V. A linear coding from transmitting a set of correlated
signals. IRE Trans. Inform. Theory, 2:41–46, September 1956.
[117] Christopher C. Paige and Michael A. Saunders. LSQR: an algorithm for sparse linear
equations and sparse least squares. ACM Trans. Math. Software, 8(1):43–71, 1982.
[119] A. G. Ramm and A. I. Katsevich. Radon Transform and Local Tomography. CRC
Press, 1996.
[121] K. R. Rao and P. Yip. Discrete Cosine Transform: Algorithms, Advantages, Appli-
cations. Academic Press, 1990.
[124] S. Sardy, A. Bruce, and P. Tseng. Block coordinate relaxation methods for non-
parametric signal denoising with wavelet dictionaries. Web page, October 1998.
[Link]
[125] S. Sardy, A. G. Bruce, and P. Tseng. Block coordinate relaxation methods for non-
parametric signal de-noising. Received from E. Candès, October 1998.
[129] Frank Spitzer. Markov random fields and Gibbs ensembles. American Mathematical
Monthly, 78:142–54, February 1971.
[130] A.S. Stern, D.L Donoho, and Hoch J.C. Iterative thresholding and minimum 1 -norm
reconstruction. based on personal communication, 1996 or later.
[131] Mann Steve and Simon Haykin. Adaptive “chirplet” transform: an adaptive general-
ization of the wavelet transform. Optical Engineering, 31(6):1243–56, June 1992.
[132] Gilbert Strang. Wavelets and Filter Banks. Wellesley-Cambridge Press, 1996.
[133] Robert Tibshirani. Regression shrinkage and selection via the LASSO. J. the Royal
Statistical Society, Series B, 58:267–288, 1996.
[134] Richard Tolimieri, Myoung An, and Chao Lu. Algorithms for Discrete Fourier Trans-
form and Convolution. Springer, 2nd edition, 1997.
[135] Paul Tseng. Dual coordinate ascent methods for non-strictly convex minimization.
Mathematical Programming, 59:231–47, 1993.
[136] Charles Van Loan. Computational Frameworks for the Fast Fourier Transform. SIAM,
1992.
[138] S. Vembu, S. Verdú, and Y. Steinberg. The source-channel separation theorem revis-
ited. IEEE Trans. on Inform. Theory, 41(1):44–54, Jan. 1995.
[139] Z. Wang and B.R. Hunt. Comparative performance of two different versions of the
discrete cosine transform. IEEE Transactions on Acoustics, Speech and Signal Pro-
cessing, ASSP-32(2):450–3, 1984.
[140] Zhongde Wang. Reconsideration of “a fast computational algorithm for the discrete
cosine transform”. IEEE Transactions on Communications, Com-31(1):121–3, Jan.
1983.
182 BIBLIOGRAPHY
[141] Mladen Victor Wickerhauser. Smooth localized orthonormal bases. Comptes Rendus
de l’Académie des Sciences de Paris, 316:423–7, 1993.
[142] Mladen Victor Wickerhauser. Adapted wavelet analysis from theory to software.
Wellesley, 1994.
[145] Song Chun Zhu, Yingnian Wu, and David Mumford. Filters, random fields and max-
imum entropy (frame): towards a unified theory for texture modeling. International
J. Computer Vision, 27(2):107–26, 1998.