
Citation 
 Permanent Link:
 http://ufdc.ufl.edu/AA00034601/00001
Material Information
 Title:
 Matrix decomposition in minimax algebra and applications in image processing
 Creator:
 Sussner, Peter
 Publication Date:
 1996
 Language:
 English
 Physical Description:
 vi, 120 leaves : ill. ; 29 cm.
Subjects
 Subjects / Keywords:
 Algebra ( jstor )
Algorithms ( jstor ) Approximation ( jstor ) Image processing ( jstor ) Integers ( jstor ) Linear algebra ( jstor ) Mathematical vectors ( jstor ) Mathematics ( jstor ) Matrices ( jstor ) Minimax ( jstor ) Dissertations, Academic  Mathematics  UF Mathematics thesis, Ph. D
 Genre:
 bibliography ( marcgt )
nonfiction ( marcgt )
Notes
 Thesis:
 Thesis (Ph. D.)University of Florida, 1996.
 Bibliography:
 Includes bibliographical references (leaves 115119).
 General Note:
 Typescript.
 General Note:
 Vita.
 Statement of Responsibility:
 by Peter Sussner.
Record Information
 Source Institution:
 University of Florida
 Holding Location:
 University of Florida
 Rights Management:
 The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. Â§107) for nonprofit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
 Resource Identifier:
 023837122 ( ALEPH )
35793543 ( OCLC )

Downloads 
This item has the following downloads:

Full Text 
MATRIX DECOMPOSITION IN MINIMAX ALGEBRA
AND APPLICATIONS IN IMAGE PROCESSING
By
PETER SUSSNER
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1996
UNIVERSITY OF FLORIDA UIBRAR!ES
ACKNOWLEDGMENTS
First and foremost, I would like to thank Dr. Gerhard X. Ritter for providing me with excellent opportunities to perform interdisciplinary research on his contract and for guiding my learning process as a mathematician. I am grateful to Dr. Panos M. Pardalos for giving new direction and inspiration to my research. I thank all my committee members, in particular Dr. David C. Wilson, for many stimulating and helpful conversations. I acknowledge Dr. Hongchi Shi, Dr. Joseph Wilson, Mrs. Monica Sweat, and all the faculty and graduate students at the Departemt of CISE who assisted me in many technical and computer science related questions. My father, Walter Sussner, deserves my deepest gratitude for his neverending support and encouragement. I am also indebted to Dr. Patrick Coffield and Dr. Sam Lambert from the Air Force Armament Laboratory for funding this research under US Air Force Contract F0863589C0134.
Finally, I acknowledge the American and the German taxpayer for their investment in my education. I hope that I will be able to justify their financial support through my future contributions to society.
ii
TABLE OF CONTENTS
ACKNOWLEDGMENTS .................................. ii
ABSTRACT ................................................ v
CHAPTERS
1 INTRODUCTION .................................... I
1.1 M otivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Introduction to Template Decomposition ..................... 2
1.3 Summary of Results and Organization ...................... 6
1.4 N otations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2 IMAGE ALGEBRA BACKGROUND ............................ 9
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Operands ............................................ 10
2.3 Operations on Images ............................... 12
2.4 Operations between Images and Templates ..................... 13
2.5 Operations between Templates .............................. 15
2.6 Metrics on Matrices and Invariant Templates with Finite Support ...... 17 2.7 Decompositions of Templates ............................... 18
2.8 Relevance of Template Decompositions ..................... 18
3 MINIMAX ALGEBRA CONCEPTS AND EXTENSIONS .............. 20
3.1 Introduction .......................................... 20
3.2 Belts and Blogs ........................................ 21
3.3 Spaces of Vectors and Matrices .......................... 26
3.4 Conjugacy ........................................... 30
3.5 The Principles of Opening and Closing ..................... 32
3.6 Linear Independence and Rank . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4 RANKBASED DECOMPOSITIONS OF MORPHOLOGICAL TEMPLATES . 45
4.1 Introduction .......................................... 45
4.2 Correspondence between Matrices and Rectangular Templates ........ 46 4.3 Invariance Properties of the Rank of Realvalued Matrices ........... 49
4.4 Separable Matrices and Templates . . . . . . . . . . . . . . . . . . . . . . . . 51
4.5 Matrices and Templates having Rank 2 . . . . . . . . . . . . . . . . . . . . . 56
4.6 Matrices and Templates of Arbitrary Rank . . . . . . . . . . . . . . . . . . . 64
iii
5 THE SVD AND THE MED .................................. 80
5.1 Introduction .......................................... 80
5.2 The Singular Value Decomposition ........................ 81
5.3 The Minimax Eigenvector Decomposition (MED) ............... 83
5.4 Comparison of the MED and the SVD ......................... 86
6 MORPHOLOGICAL DECOMPOSITION AND APPROXIMATION ...... 93 6.1 Introduction .................... 93
6.2 Morphological Decomposition and Approximation Problems ......... .93
6.3 Integer Programming Formulations ........................ 97
7 COMPLEXITY OF MORPHOLOGICAL DECOMPOSITION PROBLEMS . 105 7.1 Introduction .. . . .. . .. . . . . . . .. .. . . . .. . . . . . . .. . . . . 105
7.2 An NPComplete Problem ................................. 105
7.3 Extensions .......................................... 110
8 CONCLUSION AND SUGGESTIONS FOR FURTHER RESEARCH .... 112
REFERENCES ........................................ 115
BIOGRAPHICAL SKETCH ................................ 120
iv
Abstract of Dissertation Presented to the Graduate School of the University of Florida
in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy
MATRIX DECOMPOSITION IN MINIMAX ALGEBRA
AND APPLICATIONS IN IMAGE PROCESSING By
Peter Sussner
August 1996
Chairman: Dr. Gerhard X. Ritter
Major Department: Mathematics
Minimax algebra is a mathematical theory which originated from problems in operations research and machine scheduling. Surprising parallelisms exist between linear algebra and minimax algebra. Both linear algebra and minimax algebra can be embedded into the mathematical theory of image algebra which forms a unified mathematical environment for all image processing and computer vision applications.
Methods for matrix decomposition in linear algebra have found numerous applications in image processing, in particular for the problem of template decomposition. Therefore, it seems reasonable to investigate matrix decomposition techniques in minimax algebra with applications in image processing. We provide a mathematical basis for these investigations by establishing a new theory of rank within minimax algebra. Thus far,
V
only minimax decompositions of rank I matrices into outer product expansions are known to the image processing community. This dissertation derives various new methods for the decomposition of matrices of rank larger than 1 within minimax algebra.
Moreover, this thesis proves the NPcompleteness of the problem of decomposing arbitrary matrices in minimax algebra. Consequently, we obtain a fundamental result for the area of morphological image processing: the NPcompleteness of the general morphological template decomposition problem.
Vi
CHAPTER 1
INTRODUCTION
1.1 Motivation
A number of researchers have discovered that a variety of problems in operations research and mathematical economy, e.g. some problems in optimization on graphs and networks, in machinescheduling, and approximation theory, can be formulated using a certain algebraic lattice structure [3, 5, 8, 13, 27]. We can think of this algebraic structure as the real numbers (perhaps extended by the symbols +oo and/or co) together with operations maximum and/or minimum. In the problem formulations, which we referred to above, similarities to concepts from linear algebra become apparent.
In analogy to linear algebra, CuninghameGreen presented the minimax algebra, a unified account of the algebra of spaces of vectors over a lattice [14]. Jennifer Davidson showed how to embed the minimax algebra into the image algebra [15, 16]. The theory of image algebra incorporates all known image processing techniques in a translucent manner [19, 21, 49, 47, 52]. One powerful implication of this is that all the tools of minimax algebra are directly applicable to solving problems in image processing whenever the minimax product of matrices, denoted by [0, is involved. Jennifer Davidson's embedding of the minimax algebra into the image algebra maps vectors to images and it maps matrices to templates. This dissertation makes use of two different mappings. One of these mappings relates matrices to images, and the other one relates matrices to rectangular templates. One particular consequence of this relationship between matrices and images is that the Singular Value Decomposition known from
I
2
numerical linear algebra can be used for the purpose of image restoration and image compression [1, 2, 35]. On the other hand, matrix decomposition methods lead directly to methods for template decomposition. The problem of template decomposition  both in the linear and in the nonlinear domain  has been extensively investigated by a host of researchers. We present a brief review of different approaches to this problem below.
A multitude of researchers have applied matrix decomposition techniques stemming from linear algebra to image processing. However, the minimax algebra domain has been widely neglected in this respect. Since operations derived from minimax algebra also play an important role in image processing, we consider it timely to present an exposition on the theory of matrix decomposition in minimax algebra with applications in image processing.
1.2 Introduction to Template Decomposition
The theory of mathematical morphology grew out of the belief that many image processing operations can be expressed in terms of a few elementary operations [32, 33, 41, 54, 55]. The origins of mathematical morphology lie in the work done by H. Minkowski and H. Hadwiger on geometric measure theory and integral geometry [31, 42, 43]. Both linear convolution and morphological methods are widely used in image processing. One of the common characteristics among them is that they both require applying a template to a given image, pixel by pixel, to yield a new image. In the case of convolution, the template is usually called convolution window or mask; while in mathematical morphology, it is referred to as structuring element. Templates used in realizing linear convolutions are often referred to as linear templates. Templates can vary greatly in their weights, sizes and shapes, depending on the specific applications.
3
Intuitively, the problem of template decomposition is that given a template t, find a sequence of smaller templates t,... , t, such that applying t to an image is equivalent to applying t ,..., t, sequentially to the image. In other words, t can be algebraically expressed in terms of t1,. . . ,t'.
One purpose of template decomposition is to fit the support of the template (i.e. the convolution kernel) optimally into an existing machine constrained by its hardware configuration. For example, ERIM's CytoComputer [56] cannot deal with templates of size larger than 3 x 3 on each pipeline stage. Thus, a large template, intended for image processing on a CytoComputer, has to be decomposed into a sequence of 3 x 3 or smaller templates.
A more important motivation for template decomposition is to speed up template operations. For large convolution masks, the computation cost resulting from implementation can be prohibitive. However, in many instances, this cost can be significantly reduced by decomposing the masks or templates into a sequence of smaller templates. For instance, the linear convolution of an image with a grayvalued it x n template requires n2 multiplications and n2 1 additions to compute a new image pixel value; while the same convolution computed with an 1 x n row template followed by an n x 1 column template takes only 2n multiplications and 2(n  1) additions for each new image pixel value. This cost saving may still hold for parallel architectures such as mesh connected array processors [36], where the cost is proportional to the size of the template.
The problem of decomposing morphological templates has been investigated by a host of researchers. Zhuang and Haralick [62] gave a heuristic algorithm based on tree search that can find an optimal twopoint decomposition of a morphological template if such a decomposition exits. A twopoint decomposition consists of a sequence of templates each consisting of at most two points. A twopoint decomposition may be best suited
4
for parallel architectures with a limited number of local connections since each twopoint template can be applied to an entire image in a multiplyshiftaccumulate cycle [36]. Xu [60] has developed an algorithm, using chain code information, for the decomposition of convex morphological templates for twopoint system configurations. Again using chaincode information, Park and Chin [46] provide an optimal decomposition of convex morphological templates for 4connected meshes. However, all the above decomposition methods work only on binary morphological templates and do not extend to grayscale morphological templates.
A very successful general theory for the decomposition of templates, in both the linear and morphological domain, evolved from the theory of image algebra [19, 21, 47, 52, 48] which provides an algebraic foundation for image processing and computer vision tasks [49]. In this setting, Ritter and Gader [21, 50] presented efficient methods for decomposing discrete Fourier transform templates. Zhu and Ritter [61] employ the general matrix product to provide novel computational methods for computing the fast Fourier transform, the fast Walsh transform, the generalized fast Walsh transform, as well as a fast wavelet transform.
In image algebra, template decomposition problems, for both linear and morphological template operations, can be reformulated in terms of corresponding matrix or polynomial factorization. Manseur and Wilson [39] used matrix as well as polynomial factorization techniques to decompose twodimensional linear templates of size m x n into sums and products of 3 x 3 templates. Li [37] was the first to investigate polynomial factorization methods for morphological templates. He provides a uniform representation of morphological templates in terms of polynomials, thus reducing the problem of decomposing a morphological template to the problem of factoring the corresponding polynomials. His approach provides for the decomposition of onedimensional morpho
5
logical templates into factors of twopoint templates. Crosby [11] extends Li's method to twodimensional morphological templates.
Davidson [18] proved that any morphological template has a weak local decomposition for meshconnected array processors. Davidson's existence theorem provides a theoretical foundation for morphological template decomposition, yet the algorithm conceived in its constructive proof is not very efficient. Takriti and Gader formulate the general problem of morphological template decomposition as optimization problems [23, 58]. Sussner, Pardalos and Ritter [57] use a similar approach to solve the even more general problem of morphological template approximation. However, since these problems are inherently NPcomplete, researchers try to exploit the special structure of certain morphological templates in order to find decomposition algorithms. For example, Li and Ritter [38] provide very simple matrix techniques for decomposing binary as well as grayscale linear and morphological convex templates. A separable template is a template that can be expressed in terms of two onedimensional templates consisting of a row and a column template. Gader [22] uses matrix methods for decomposing any grayscale morphological template into a sum of a separable template and a totally nonseparable template. If the original template is separable, then Gader's decomposition yields a separable decomposition. If the original template is not separable, then his method yields the closest separable template to the original in the mean square sense.
Separable templates are particularly easy to decompose and the decomposition of separable templates into a product of vertical and horizontal strip templates can be used as a first step for the decomposition into a form which matches the neighborhood configuration of a particular parallel architecture. In the linear case, separable templates are also called rank 1 templates since their corresponding matrices are rank 1 matrices. O'Leary [44] showed that any linear template of rank 1 can be factored exactly into a
6
product of 3 x 3 linear templates. Templates of higher rank are usually not as efficiently decomposable. However, the rank of a template determines upper bounds of worstcase scenarios. For example, a linear template of rank 2 always decomposes into a sum of two separable templates.
1.3 Summary of Results and Organization
In nonlinear or morphological image processing, convolutions of an image with a template consist of combinations of maximum, minimum, and addition operations. The theory of minimax algebra provides the mathematical background for the algebraic structures arising from these operations. This dissertation presents a new concept of matrix rank within minimax algebra which is similar to the one in linear algebra. Due to the correspondence of matrices and rectangular templates, we obtain a concept of morphological template rank. In this context, separable morphological templates represent exactly the templates of rank 1. A separable template is easy to decompose into an outer product of a column and a row template. The cost savings in computational complexity achieved by a decomposition of this form still remain, if templates of arbitrary, but fixed rank k are decomposed into outer products of column and row templates. First, we present an algorithm which, for any rectangular morphological template of rank 2, constructs a weak decomposition of this template into two column templates and two row templates in polynomial time. Afterwards, we present a heuristic algorithm for the decomposition of arbitrary templates of relatively small rank k.
In linear algebra, an exact decomposition of a rank k matrix into k matrices of rank 1 can be determined by using the SVD. This thesis provides a matrix decomposition in minimax algebra which is mathematically very similar to the SVD in linear algebra. However, this decomposition which we call MED, generally does not provide the best
7
approximation of a matrix by k outer products of column and row templates in terms of an appropriate matrix norm. Instead, we solve this problem by formulating it as an integer programming problem. Furthermore, the general morphological template approximation problem can be formulated as an integer programming problem. We provide preliminary computational results.
Finally, we justify these integer programming approaches by showing that the general morphological template decomposition problem is NPcomplete. Indeed, the problem of weakly decomposing arbitrary rectangular templates into column and row templates is already NPcomplete.
This dissertation is organized as follows: First, we introduce the reader to the language of image algebra. This mathematical theory is suited to describe all image processing operands and operations in a translucent manner. Since we focus on morphological image processing, we proceed with a brief description of the algebraic structures defined in minimax algebra. We relate matrices in minimax algebra to rectangular morphological templates. In chapter 3, we establish a rank method for the morphological decomposition of matrices and rectangular templates. We briefly review the Singular Value Decomposition in chapter 4. We derive a similar decomposition in minimax algebra, called Minimax Eigenvector Decomposition. In chapter 5, new integer programing approaches for solving arbitrary instances of the general morphological template decomposition problem are presented. Finally, we justify this strategy by showing that the problem of decomposing arbitrary morphological templates is NPcomplete.
1.4 Notations
Before we proceed, let us discuss some frequently used notations. As usual, the symbol N stands for the set of natural numbers, Z for the set of integers, and R stands
8
for the set of real numbers. The symbol R+ denotes the positive real numbers. The real numbers extended by the symbol oo are denoted by R, and the real numbers extended by the symbol +oo are denoted by R+,. Similarly, R+. stands for R+ U {oo}, and R+O stands for R+ U {+oo}. When adjoining both +oo and oo to R, one obtains the set R ,,.
We use the symbols V and A to denote the binary operations of maximum and minimum. The positive part of a real number x, denoted by x+, is defined as 0 V x. Similarly the negative part of x is x = 0 V (x).
Matrix inequalities are always defined entrywise. For example, given matrices A and B, the inequality A < B means ai3 < bij V i, j.
If y is an associative and commutative operation on a finite set X = {x1,... x,}, then we abbreviate the expression x1... 7X by IF 1(xi) or by F xi. The symbol
i= iE1t ,.,n}
F x is yet another way to denote x1Y... 7X provided that x denotes the transpose of (x1,.. .,xn), or briefly (xi,.. .,xn)T. If f is a function on X, then f(X) denotes {f(x1),. . . ,f(Xn)}. A sequence (ki,k2,. ..,k) is called strictly increasing if k, < k2 < ... < kn.
CHAPTER 2
IMAGE ALGEBRA BACKGROUND
2.1 Introduction
Images not only provide important sources of information in every day life, but also in various sciences and in the military. However, often large parts of the information contained in an image remains hidden from the eye of the human observer. Image processing techniques, when implemented on a computer, can help to filter out these parts of the information. For instance, image processing methods can be applied in order to identify and locate objects in an image. The objects the user is looking for can be as diverse as cancer cells in medical applications and tanks in military applications. In cardiology, for example, physicians need to find the location of the anterior and posterior wall of the heart.
A multitude of image processing methods for image smoothing, image enhancement, edge detection, thinning, skeletonizing etc. exist [4, 30, 29]. When faced with a specific problem, like finding certain objects in an image, the developer of an image processing algorithm has to make the difficult decision which techniques to use. In this case, a certain amount of a priori knowledge, e.g. on the shape and size of the objects, will turn out to be very useful.
For a long time developers of image processing software have documented their algorithms by using lengthy word descriptions and analogies that are extremely cumbersome and often ambiguous. These ad hoc approaches led to a proliferation of nonstandard notation and increased research and development cost. In response to this chaotic situation,
9
10
G. X. Ritter et al. have developed the (AFATL) image algebra since the mid1980's. This algebraic structure is specifically designed for image manipulation.
The image algebra constitutes a heterogeneous or manyvalued algebra in the sense of Birkhoff and Lipson [7, 48]. Intuitively, an algebra is simply a collection of nonempty sets together with finite number of mappings (operations) from one of the sets into another. In image algebra, six basic types of operands exist, namely value sets, coordinate sets, the elements of each of these sets, images and templates. In this dissertation, we restrict our attention to singlevalued images. We refer the reader to other publications for a discussion of multivalued images [48, 52].
2.2 Operands
A point set is a subset of a topological space. In digital image processing, point sets are often assumed to be rectangular discrete arrays. A value set is a homogeneous algebra such as Z, R, R. (= R U {co}), R +(= R U {oo}) or C.
An Fvalued image a is the graph of a function a : X * F where X denotes some point set and F denotes some value set. Obviously, a can be represented as a matrix if X is a rectangular discrete array. The set of all images a X * F is denoted by FX.
x
The notation RX0 is frequently used instead of (R,)
Let X, Y be point sets and F a value set with neutral element 0. A generalized Fvalued template t from Y to X is a function t : Y > FX. The set of all templates t : Y > FX is denoted by (FX)y. So for each y E Y the value t(y) is an Fvalued image on X. We define ty = t(y) and call the point y the target point of t. The function values ty(x), where x E X, are called weights of the template t at location y.
I1
The support of ty is defined as follows:
S(ty)= {x E X : ty(x) }
Note that the neutral element 0 depends on the algebraic structure of F. For example, the element 0 represents the neutral element for R together with +, whereas the element
oo represents the neutral element for R, together with V.
A template t E (F) X is called translation invariant if and only if for each triple x, y, z E X with x + z C X and y + z E X the following equation holds:
ty(x) = ty+z(x + z) .
In this case, it suffices to know tyO for an arbitrary yo E X. If X is discrete and S(tyo) is finite then t can be represented pictorially.
For example, let X = Z2. Moreover, let r E (R3,)X be the translation invariant template which is determined at each point y = (x, y) C X by the following function values of x E X:
7 if x=y
4 if x= y (0, 1)
ry(x)= 12 if x = y +(0, 1)
2+ry(x+(1,0)) if xi =x1 and y1
o else .
We can visualize the invariant template r as shown in Figure 1.
12
y1 y y+1
X1 2 5 10
x ry 4 7 12
3 0 5
Figure 1. The support of the template r at point y. The hashed
cell indicates the location of the target point y = (x,y).
A translation invariant template r is called rectangular, if S(ry) forms a rectangular discrete array. We speak in particular of an m x n template, if S(ry) has size m x n. An m x 1 template is also called a column template of length m, and a 1 x n template is called a row template of length n. For the remainder of this dissertation, point sets can always be identified by the symbols X, Y, Z and value sets by the symbols F, G, H. The point set X is always assumed to be finite.
2.3 Operations on Images
The basic operations on Fvalued images stem from the naturally induced operations on the value set F. In particular, addition, multiplication, and maximum are defined as follows:
a + b = {(x, c(x)) : c(x) = a(x) + b(x), Vx E X} V a, b E R,
a  b ={(x, c(x)) : c(x) =a(x) . b(x), Vx E X} V a, b E Rx
a V b {(x, c(x)) : c(x) = a(x)Vb(x), Vx E X} V a, b E Rx, . We adopt the usual convention r + (oo) = (oo) + r = oo for all r E R 0.
13
2.4 Operations between Images and Templates
The more complex operations which we define in this section are based on a binary operation o : F x 6 + H and an associative and commutative operation 7 on the value set H.
If a E GX and t C (FY) X, then a 6 9t E HY  the forwards generalized product of the image a with the template t  can be formed:
a~t= {(y, b(y)) : b(y) =LF (a(x) o tx(y)), y E Y
Hence, a forwards product of this form can effect a transformation from one image having a certain value set and a certain point to another image having a different value set and a different point set. This thesis concentrates on the backwards products of images and templates which are defined below.
If a E FX and t E (GX) Y, then a Ot E HY  the (backwards) generalized product of the image a with the template t  is defined as the following image:
a Ot = (y, b(y)) : b(y) =F (a(x) o ty(x)), y Y}.
Specific examples of this operation are given below. Unless stated otherwise, we are referring to backwards imagetemplate products for the remainder of this thesis. If, for example, a E Rx0 and t E (R0) Y, then the generalized convolution ("ï¿½") and the additive maximum (" M ") are defined as follows:
a M t = (y, b(y)) : b(y) = Y a(x) + ty (x), y C Y xEX )}
" (Dt = (y, b(y)) :b(y) = 1:a(x) * ty(x), y E Y xEX
14
where
V a(x)+t,(x)=maxia(x)+t,(x),xcX. xEX
The operation additive minimum (" E ") is defined in a similar fashion.
Many image processing techniques such as the Rolling Ball Algorithm and algorithms for noise removal employ imagetemplate products of the form a M t or a 01 t. Both forms appear in the MaxMin Sharpening Algorithm, which enhances blurred images.
Figure 2. Blurred characters and result of applying maxmin sharpening to blurred characters
One of the most widely used transforms in image processing is the Fourier transform [10, 59]. Like many other image processing transformations, the Fourier transform can be expressed in terms of a generalized convolution. This means that there is a template f such that the Fourier transform of an image a is given by a Gf.
The images in Figure 4 show an image of a jet and its corresponding Fourier transform image after centering. Notice that almost all of the information is contained in the center of the latter image.
AeJ_
15
I
Figure 3. The image of a jet and its centered Fourier transform image.
Typical applications of the Fourier transform are low and high pass filtering [29, 49]. The figures below demonstrate how low pass filtering can be applied to medical images.
Figure 4. Noisy angiogram and filtered angiogram using a low pass filter
2.5 Operations between Templates
The basic operations between generalized templates are induced by the basic operations on images. Thus, the addition, the multiplication, and the maximum of two templates are straightforward:
s + t ={(y, ry) :ry =sy +ty, y E Y} V s,t E (Rx)
s * t ={(y, ry) :ry =sy * ty, y E Y} V s,t E (RX
S V t = (y, ry) :ry =Y syv ty, y E Y} V s, t E (RXO)
16
Since the set R_, together with the morphological operations + and V forms an algebraic structure in minimax algebra, the R,valued templates are also called morphological templates.
Templates can not only be convolved with images, but also with other templates. Ritter [48, 52] defines the generalized product of a template s E (FZ) X and a template t C (GX) Y, also denoted by s Ot, as the template r E (HZ) Y by setting
ry(z) =.F (sX(z) a ty(x)) Vy C Y, Vz E Z.
For example, realvalued templates t and s allow us to form a generalized convolution (" ï¿½") and an additive maximum (" M "):
(s ï¿½t) (z) = 1: (sX (z) * ty (x)) Vy E Y, Vz E Z xEX
(s M t)y(z) = V (sx(z) +ty(x)) Vy E Y, Vz E Z xEX
The following column templates r, s, t C (R00) X satisfy r s Mt.
2
3
r 4 Sy= 0 Y 2
5 3 3
6
Figure 5. The template r constitutes the additive maximum of the templates s and the template t.
Note that the generalized products s Ot and t Os of an m x 1 rectangular template s and a 1 x n rectangular template t always lead to m x n rectangular templates.
17
Suppose that X = Z2 and that s, t E (Rx)X are the following translation invariant templates: The template product r = s M t is given by the template r in Figure 1.
2
SY 0 t= 7
7
Figure 6. Pictorial representation of a column template s and a row template t.
2.6 Metrics on Matrices and Invariant Templates with Finite Support
It is well known that a metric on a metric space M induces a metric on the space M1" of all m x n matrices. Let d be the metric on R given by: d(a, b) =ja  bj Va, b E R .
Then the following definitions induce metrics on Rmxn: n M
di(A, B) = d(aij, bij)
=1 3=1
d2(A, B) = V d(a i, bij) .
2,)j
Similarly, if s E (RX)Y and t E (Rx) are invariant templates with the same finite support S at some target point y E Y, then the following equations provide metrics on the set of all invariant templates with that particular finite support: d,(s, t) = d(ty(x), sy(x)) xES
d2(s, t) = y d(ty(x), sy(x)) .
xES
Since the templates s and t are invariant, these definitions do not depend on the choice of the target point y E Y.
18
2.7 Decompositions of Templates
A sequence of templates (t,....., t'") in (FX) X together with a strictly increasing sequence of natural numbers (ki,. . . , k,) is called a weak decomposition (with respect to the operation "0") of a template t E (FX)X if the template t can be represented as follows:
t = (t (otk. (tk. +1... o)tk2>. (tkn.+l 0... oDtkn)
In the case where n = 1, we speak of a strong decomposition of the template t. For example, a strong morphological decomposition of a template t E (R)S)X has the following form:
t = t, MOt2M ..M t
If the template t strongly decomposes into a column template tV and a row template t2, it is referred to as a separable template.
For example, the template r E (R&O )X given in Figure 1 represents a separable template since this template decomposes into the column template s E (Rxo)X and the row template t E (R,)X of Figure 6.
2.8 Relevance of Template Decompositions
The following associative and distributive laws hold for an arbitrary image a E FX and arbitrary templates t E Fx and s E Fx: a M (s M t) = (aM s) M t,
a M(sVt) = (a Ms)V (a Mt) .
19
These results establish the importance of template decomposition.
Realizing the convolution of an image a with a template t on a computer is an extremely complex task, since each pixel value at location y depends on all pixel values of a within the neighborhood given by ty. Decomposing the template t often helps to reduce the number of computations involved in forming an imagetemplate product by providing alternative ways to compute this convolution. For example, the application of a separable template to an image is equivalent to the sequential applications of a column template and a row template. Hence, a quadratic problem reduces to a linear problem. Computing the Fourier Transform of an image represents a linear convolution. The speedup of this computation achieved by the Fast Fourier Transform relies heavily on linear template decomposition.
Template decompositions have also been successfully used for the purpose of achieving a better fit to the computer architecture. For example, if each processing element has direct neighbors in N, S, W, and E directions, then the user of image processing software benefits from template decompositions into invariant templates having supports of the following form:
77
CHAPTER 3
MINIMAX ALGEBRA CONCEPTS AND EXTENSIONS
3.1 Introduction
In analogy to linear algebra, CuninghameGreen presented the minimax algebra, a unified theory of the algebra of spaces of vectors over F [14]. Minimax algebra includes familiar concepts from linear algebra such as linear dependence of vectors, matrix rank, and eigenvalues of matrices. Operations on vectors and matrices in minimax algebra resemble the operations on vectors and matrices in linear algebra. In his Ph.D. dissertation, Paul Gader embedded the linear algebra into the image algebra [20]. Likewise, J. L. Davidson provided an isomorphism of the minimax algebra into the image algebra [15, 17]. This isomorphism maps vectors to images and matrices to templates. Hence, results on matrix decomposition in minimax algebra induce results on template decomposition in image algebra.
This dissertation exploits a different kind of relationship between matrices in minimax algebra and rectangular morphological templates, which is also useful for the purpose of template decomposition. A similar relationship exists between matrices in linear algebra and rectangular linear templates [48]. In this context, the rank of a rectangular realvalued template is given by rank of the corresponding realvalued matrix. In analogy to the linear domain, this thesis defines a notion of matrix rank within minimax algebra which leads to a notion of rank of rectangular morphological templates.
20
21
3.2 Belts and Blogs
Consider an algebraic structure (B, *) consisting of a set B together with a binary operation * which assumes the role of addition in B. We call (B, *) a commutative band or semilattice if the following axioms hold for all x, y, z E B:
X1 :x*y = yOx .
X2 :x O(y OZ) =(X Oy) OZ.
X3 :xox = x
Homomorphisms between commutative bands (B1, o) and (B2, o) are considered to be functions f : B1 + B2, which satisfy the relation f(x o y) = f(x) c f(y). Although denoted by the same symbol o, the addition operations in B1 and in B2 are not necessarily the same.
Usually, the operations o can be thought of as either the maximum operation or the minimum operation. For these interpretations of o, an additional axiom Xo holds for all x, y E B:
Xo : x o y = x or xoy = y.
The axioms XI, X2, and X3 induce a partial order on B by means of the following definition:
x>y. x =x y Vx,yGB.
Recall that a set B is partially ordered if and only if the following axioms hold for all x, y, z E B:
Y1: x > y and y > x > x = y
Y2: x > y and y > z > x > z
22
Another partial order on a commutative band B is given by the relation <, where x < y if and only if y > x. We write x > y or y < x whenever x > y and x 5 y.
A commutative band (B, *) satisfying the additional axiom Xo is called linear. In this case, the relation > establishes a linear order on B and satisfies an axiom Y corresponding to the axiom Xo.
YO : x > y or y > x .
The real numbers R as well as the extended real numbers R .. together with V or A are examples of linear commutative bands.
A belt is an algebraic structure (E, o, *) consisting of a set E and two binary operations o and *. The operation o represents the addition in E, and the operation
* represents the multiplication in E. The system (E, o, *) satisfies the following axioms (x, y, z C E):
X1, X2, X3 : (E, o) is a commutative band.
X4 : x*(y*z)= (x*y)*z.
X5 x*(yOz) = (x*y)cO(x*z) X6 (yOz)*x = (y*x)cO(z*x)
Note that the algebraic structure (E, o, *) is known to be a semilatticeordered semigroup. However, since this term seemed so unhandy, CuninghameGreen coined the term belt, evoking a little flavor of the words ring and band.
Like in linear algebra, there is a natural way to define homomorphisms between algebraic structures such as belts. Given two belts (E1, 0, *) and (E2, o, *), a belt homomorphism from E1 to E2 is a function f : E1 * E2 which satisfies the conditions stated below. The operations of addition and multiplication in E1 are not to be confused
23
with addition and multiplication in E2.
f(xcOy) = f(X) Of(y),
f(x * y) = f(x) * f(y)
Vx,y E E1 .
If the commutative band (E, o) underlying a belt (E, o, *) is linear, i.e. if this algebraic structure also satisfies axiom X0, then the belt (E, o, *) is called linear as well. We speak of a commutative belt whenever the operation * is commutative, i.e. whenever the following axiom X7 holds:
X7: x * y = y * x V x,y E E.
A belt satisfying additional axioms X8 and X9 will be called a division belt.
X8 : 0 E E such that 0 * x = x *0 = x V x E E.
X9 : V x E E 3 x' E E such that x' *x = x * x' = 0.
Axioms X8 and Xg establish the existence of an identity element and of inverse elements with respect to the multiplication *. As usual, the multiplicative inverse of x C E is denoted by x1. The real numbers together with V and +, denoted by (R, V, +), form a division belt. It is easy to recognize that the following systems belong to the class of linear belts, but not to the class of division belts: i. (R00, V,+)
ii. (R+,O, V,+)
iii. (R+,,, V,7)
iv. (R+ m, V,+)
24
Lemma 3.2.1. The following properties hold for arbitrary elements x, y, z of a linear division belt (E, o, *).
1.
X > Y X 0 z > y 0 z x*z >y*z and z*x > z*y.
2.
x > y => x*z > y*z and z*x > z*y.
3.
z>x and z>y > z>xoy.
Proof. For a proof of Property 1, we refer to Proposition 22 and Proposition 25 in CuninghameGreen's Minimax Algebra.
Suppose that x > y for some x, y E E. By Property 1, the inequality x * z > y * z holds for all z E E. The assumption x * z = y * z leads to an immediate contradiction, if we multiply by z1 on the right. Hence, x * z > y * z as desired. Similary, we obtain z * x > z * y.
Let z > x and z > y for some x, y, z E E. Consider z o (x o y). By the axiom of associativity X0, z o (x o y) = (z o x) * y. Using the linearity of the belt (E, 0, *) and the facts that z > x and z > y, we obtain that z o (x o y) = z, but z o (x o y) 5 x o y. In other words, z > (X 0 y).
25
Blogs are progressively constructed from division belts (E, o, *). First a dual addition o' is introduced for all x, y E E.
x *' y = (x )1.
Next we define a set 0 by adjoining universal bounds oo and +oo to E, i.e. these new elements satisfy
(c0) 0 x = x c (c0) = x (+0o) 0 x = x C (+o)= + (cc) >' x = X c' (cc) = cc (+co) 0' x = x 0' (+00) = x VxCO.
Then we extend the domain of the operation * to 0 x 0 as follows: (00) * (+c) = (+c) * (cc) = cc (oo) *x =x *(oo)= oo V x c E U {oo} (+oo) *x =x *(+oo) +oo V x E E U {+oo}.
Finally, we introduce a dual multiplication *' which acts exactly like * except that we stipulate
(00) *' (+cc) = (+cc) *' (cc) = +c .
Hence, a blog is given by the system (0, 0, *,0 ', *'). Note that the systems (0, o, *) and (0, *', *') form belts. The set E is called the set of finite elements of this blog. If (E, o, *) forms a linear division belt, then we speak of a linear blog. Similarly, (E, C, *) forms a commutative division belt, then we speak of a commutative blog.
26
By a direct consideration of cases, we can confirm that the following "associative inequalities" hold in an arbitrary blog (0, c, *,c", *').
Lemma 3.2.2. If (0, o, *,', *') is a blog, then all x, y, z E 0 satisfy x * (y *' z) 5 (x * y) *' z ,
x *' (y * z) 5 (x *' y) * z
The principal interpretation is (Roo, V, A, +, + . The operations + and +' act like the usual addition with the following exceptions.
(o) + (+00) (+00) + (00) = 00 (00) +' (+0) = (+00) +' (00) = +00 . CuninghameGreen shows that the only finite blog is the threeelement blog representing the set {oo, 0, +oo} together with the operations V, A, +, +' restricted to this set.
3.3 Spaces of Vectors and Matrices
In analogy to the concept of vector space over a field in linear algebra, CuninghameGreen defines a concept of bandspace, or briefly space, over a belt. Suppose (B, c) is a commutative band and (E, *, *) is a belt. Furthermore, assume that a right multiplication of elements of B by elements of E exists, i.e. there is a binary operation *: B x E + B. Then we say (B, *) is a (right) bandspace over (E, o, *), or briefly B is a space over E, if the following axioms hold for all x, y E B and for all A, p E E:
S1 x*(A*p) = (x*A)* p.
S2 (x y) *A = (x *A) (y *A) S3 : * (Aq p) = (x *A) (x* p)
27
Additionally, if E has an identity element 0 with respect to the multiplication *, we assume for all x E B:
S4 : x * 0 = x
Clearly, we can formulate "left" variants of axioms Si to S4 and define thereby a left bandspace over E. We will call B is a 2sided space over E, if B is a right bandspace over E as well as a left bandspace over E, and if the following associative law holds for all x E B and for all A, p E E:
A * (x ) (A* x) .
For any given belt E and integers m and n, the set of m x n matrices with entries in E, which is denoted by E", forms a 2sided space over E. In this setting, the addition o in E" as well as the scalar multiplication are performed entrywise. Likewise, the set E of ntuples with entries in E represents a 2sided space over E. Addition and scalar multiplication are induced by addition and multiplication in E.
The algebraic operations o and * can also be used to define a matrix product
* : EmxP x EPX"  E x", where p E N. In analogy to the usual matrix product known from linear algebra, multiplying a matrix A E EmxP by a matrix B E EP'" results in a matrix C E Emx", whose entries are given by the following formula:
c = (aj,* bj)
(31)
Vi = 1,..., m; Vj =1,.. ., n .
In particular, if E equals the belt (R ,, V, +), the matrix C representing the matrix
28
product of A and B is computed as follows:
P
cij= (air + brj)
r=1
V i = 1, ..., m; Vj = 1, ..., n. In this case, we denote the matrix product by the symbol M . Similarly, the matrix product, which is derived from the operations A and +' in the belt (Rioo, A, +'), is denoted by the symbol E .
A matrix product u * v of a column vector u and a row vector v is called an outer product. We say that a matrix A decomposes into k vector pairs, if A can be written as a
k
sum of k outer products. For the matrix product M, we have A = V (ul M v'), where i=1
u1 are column vectors and vi are row vectors for 1 = 1,..., k.
Example.
4 ) M (3 8 7)
1 9 0 =
0 5 90
[(') M(3 8 7)) V 9 ( 9 0) =
0 5
9 18 9 .
5 14 7
Theorem (CuninghameGreen) 3.3.1. The following associative inequalities hold for matrices over R oo.
Xz1(Y MZ) ;> (XEIY )M Z, X M(Y E Z) < (XfOY )& Z
V X ER ,VY E RO ',V Z E RaoxP
29
Proof. Using distributivity in the belt (R 0, V, +) and Lemma 3.2.2, we infer for arbitrary i C {1,...,m} and j E {l,...,p}
o nl
((X 9 Y) M Z) V A (Xik +' ykh) + Zh, h=1 k=1
0 71
V A[/\ (Xk + Ykh) + Zhil h=1 k=1
0 n
iV A [Xi (Ykh + Zhi)h=1 k=1
The principle of opening implies that V A [Xik +' (Ykh + Zh,)] h=1 k=1
V [Xik +' (ykh + zhi)] for all k = 1,...,n. By applying the principle of h=1
closing and distributivity in the belt (R i,, A, +'), we obtain
0 n n 0
V A [Xik 1 (Ykh + zhj)] i AV[Xk (Ykh + ZhDI
h=1 k=1 k=1 h=1
n o
[ik 1 V (Ykh + Zhj) k=1 h=1
= (X I (Y M Z))ij.
For a proof of the second half of the theorem, we refer to the proof of Theorem 63 in CuninghameGreen's Minimax Algebra.
The matrix product defined in Eq. (31) can be viewed as a special case of the generalized matrix product [9]. In conjunction with the matrix sum 0 and the matrix product *, the set E," forms a belt. Direct verification of the relevant axioms shows that EPX" is a commutative band as well as a right bandspace over the belt E"'".
Another important class of spaces over a given belt (E, o, *) is the class of function spaces, that is to say spaces in which the commutative band (B, 0) is the commutative band (EU, o) of all functions from some given set U to E. Here the addition in EU and the scalar multiplication on EU are induced by the addition and the multiplication
30
in E. Specifically, we obtain the following definitions for all u E U, for all f,g E EU, and for all A E E:
(f g)(u) = f(u) 0 g(u),
(f * A)(u) = f(u) * A,
(A * f)(u) = A * f(u)In view of these definitions, it is easy to see that the function space EU is a 2sided space over E. For example, we can take E to be the real numbers.
Although not defined in his book Minimax Algebra, CuninghameGreen suggests that bandspace homomorphisms are functions F from one bandspace B over a given belt E to another bandspace C over E which satisfy the following criteria: F(b b2) = F(bi) * F(b2) V bi, b2 E B F(A *b) = A * F(b) V A c E, V b E B. The notions isomorphism, endomorphism, and epimorphism have the usual meaning whenever they refer to homomorphisms between commutative bands, belt homomorphisms, or bandspace homomorphisms. If M and N denote the sets {1,... , m} and { 1,... ,n} for given integers m and n, then the function space EMxN can be identified via an isomorphism with the space of m x n matrices over E.
3.4 Conjugacy
In analogy to conventional linear operator theory, the theory of minimax algebra gives rise to a theory of conjugacy. Taking an axiomatic conjugacy approach, CuninghameGreen defines a certain involution x + x* subject to certain axioms. In the present section, we discuss axiomatic conjugacy for belts and matrices over belts.
31
A commutative band B2 is said to be conjugate to a given a commutative band B1, if B1 is band isomorphic to B2. Suppose now that E1 and E2 are belts. The belt E2 is conjugate to the belt E1, if there exists a belt isomorphism g from El to E2. Note that conjugacy is an equivalence relation, in particular El is conjugate to E2 with respect to the isomorphism g1. Hence, we can say that E2 is the conjugate of El. Vice versa, E1 is the conjugate of E2. We denote the conjugate of a given algebraic structure E by the notation E*. Hence (E*)* = E.
Just as in conventional linear operator theory, we shall henceforth employ the notation x* to denote the element of E* corresponding to x E E under the isomorphism from E to E*. Thus (x*)* = x. We call the element x* C E* the conjugate of x E E. A straightforward verification yields that every division belt is conjugate to itself, or briefly selfconjugate, under the isomorphism given by x* = x1. We say that every blog (0, 0, *, c', *') is selfconjugate, because the belt (0, o', *') is the conjugate of the belt (0, o, *) under the following isomorphism: (oo)* =+00
(+oo)* = 0
x* = 1,
X =X
whenever x is a finite element of 0.
Recall that, for any integers m and n, the set of m x n matrices over a belt E forms a commutative band and the set of n x n matrices over a belt E forms a belt. Hence, it is possible to consider the conjugate systems of (Emxn, .) and (Enxn, 0, *). In fact, Theorem 7.4. of Minimax Algebra shows that the conjugate systems (E"xn)* and (Enxn)* are exactly the systems (E*)nx" and (E*)nxf.
A matrix A C E" x' is mapped to its conjugate matrix A* E (E*)nx"" under the
32
isomorphism of commutative bands from (Emxn, 0) to ((E' n)*, o). The matrix A* is obtained from A by transposing and conjugating, i.e. for all i = 1,...,m and all j = 1,..., n, the (i,j)th entry of A* is equal to the conjugate of the (j, i)th entry of A. Formally, we write:
(a*)i = (aji)*
The same relationship as above holds for the belt isomorphism from (E"xn, c, *) to ((Enxn)*, 1 , *).3.5 The Principles of Opening and Closing
The principles of opening and closing provide a set of rules which are useful for the purpose of proving theorems in minimax algebra. Let (W, 0, *, *', *') be a blog, and let the symbol V denote the set of finite elements of this blog. Recall the definition of the dual addition o' on the set of finite elements of W: xc0'y = (X1o Y1[1
This equation implies that the following equation holds for all X1, X2,... Xn E V: [on(Xj)]l = (0,) I(xj), (32)
We now discuss some of the basic formal properties of the operators 0 and 0'. We first address the principle of closing.
Suppose that X1, x2,   , Xn and Y1, Y2,   . ,yn are elements of the blog W. Let the following inequality be given: Xi;>y; VI=1,...,n
(33)
xi 0 yi = x V i=..n .
33
By using the commutativity and associativity of the operations of addition c and dual addition *', we can show that:
(34)
( n') I (Xi) (O'?1Yi)Since the operation * is idempotent, i.e. x c x = x for all x E W, it may be possible to use the following facts in order to simplify the inequalities in Eq. (34):
X; = X V i =,.,n > 1 O_(X;) = (,') (;
=ox V i = 1,...,n > K 1(yi) = (, ( )n
The manipulation principle used in order to obtain Eq. (34) from Eq. (33) is called the principle of closing. This principle is expressed in Proposition 31 of Minimax Algebra:
Theorem 3.5.1. If W is a blog, then an inequality, governed by a universal quantifier, may be replaced by an inequality having the operator K (respectively 0') governing each side of the original inequality, with the index in both cases running over the same finite subset of W as for the universal quantifier.
When applying the principle of closing, we add the operator 0 or the operator K' to both sides of an inequality.
CuninghameGreen also introduces the principle of opening, which, informally speaking, relies on the following facts:
1. An expression governed by an operator Q is decreased, if the operator is removed. 2. An expression governed by an operator K' is increased, if the operator is removed.
For example, let us consider the expressions O U1(x) and (O')L1(x), where x1, X2,. ,x E W. The associativity and commutativity of the operations o and o'
34
yields two sets of inequalities:
x0 (') (i) Vi=1,...,n
The validity of these inequalities is preserved when applying identical sequences of additions, dual additions, left and right multiplications, and left and right dual multiplications on both sides. These considerations lead to the following theorem:
Theorem 3.5.2. Let an algebraic expression 0 be obtained from a wellformed algebraic expression a by deleting an operator K which does not form part of the argument of any inversion operation ()'. Then the inequality a > 0, followed by a universal quantifier whose index runs over the same (finite) set as for the operator K, is valid in W.
Dually, let an algebraic expression / be obtained from a wellformed algebraic expression a by deleting an operator 0' which does not form part of the argument of any inversion operation (K . Then the inequality a < #, followed by a universal quantifier whose index runs over the same (finite) set as for the operator 0', is valid in W.
The application of either half of this theorem will be called opening with respect to the index in question.
Example. Let a be the algebraic expression 0, =(xj) * [1'j(xj)]I. Deleting the operator K from the first part of a results in an algebraic expression / which equals Xj * [0"j1(xi)] _. Hence, the principle of opening yields:
K =1(Xj) * [KUI1(Xi)F 1 >X * [K ((Xi) Vj = 1,..
35
3.6 Linear Independence and Rank
In linear algebra, the rank of a matrix is defined as the number of linear independent row vectors or column vectors. In this section, we present CuninghameGreen's definition of matrix rank in minimax algebra, which is based on the concept of strong linear independence. Then we provide our own definition of matrix rank and relate this notion to the one given by CuninghameGreen. As a matter of convenience, we assume that (W, c, *, o', *') is a commutative and linear blog. We denote the latticeordered group of the finite elements of W by the symbol V.
Definition. A vector v E W' is said to be linearly dependent on the vectors V, .. . kv E W' if and only if there exist ci C W for all i = 1, ...,k such that
= = (ci * vS)
We define linear independence as the mere negation of linear dependence. Vectors V1, ... , vk E W" are linearly independent if no one of them is linearly dependent on the others.
Example. Let the following elements of RV be given: v 5 5 1
v= 1 ,V= 2 ,2
7 2 92
Since v = [2 + vi] V [(2) + v2], the vector v is linearly dependent on v' and v2
The book Minimax Algebra by CuninghameGreen provides a convenient mechanical procedure, called the %test, in order to check if given vectors vi, ... , vk E W" are linearly independent. The 2test involves the formation of a matrix 2 E W n"". If ai3
36
denotes the (ij)th entry of %, and if V E Wnxk is the matrix, whose jth column is vi for all j = 1,..., k, then the matrix Q is defined as follows:
{ 0 else
The 21test compares each column of V with the corresponding column of V * [, and comes to the following conclusion:
Theorem 3.6.1. For each j = 1,... , k, the jth column of V * is identical with vi if and only if vi is linearly dependent on the other columns of V E Wnx k. The elements of the jth column of % then give suitable coefficients to express the linear dependence.
Example. We illustrate this technique with an example from the principal interpretation R o, V, A, +, +). Let v", ... ,v4 be the following vectors in R3
)1=, V2 ()V3 ()V4
The next step consists of forming the matrix V* 1 V, with vi, ... ,v4 representing the column vectors of V:
1 3 o 1 3 2 3)
V*EVg = 3 4 5 E 3 4 2 1
2 2 5 o 5 5 3
3 1 3)
oo 1 1 2
oo 0 2 3
oo 0 0 2
00 0  1 0
37
When overwriting the main diagonal of this matrix with oo, we obtain the matrix 2t E Rx. We proceed by forming the product V M 2.
00 1 1 2
o o 2 3
00 0 oo 2 '
00 0 1 00)
oo 4 2 0
VM = (oo 4 2 0).
00 5 3 3
A comparison of V with V M % shows that v2 is linearly dependent on v', v3, and v4 with coefficients occurring in the second column of %, i.e.:
V2=(1+vI)V(0+v3)V(0+v4).
In conventional linear algebra, the concept of linear dependence provides a basis for a theory of rank and dimension. In minimax algebra, however, some unexpected anomalies occur regarding the notion of linear dependence. CuninghameGreen addresses this situation in the following two theorems (Theorems 164 and 165 of Minimax Algebra).
Theorem 3.6.2. Suppose W is a blog other than ({oo, 0, oo}, V, +, A, +'). If n > 2 and k > 1, then there exist k finite ntuples, no one of which is linearly dependent on the others.
Theorem 3.6.3. If W is the threeelement blog ({oo, 0, oo}, V, +, A, +'), then there exist (n2  n) finite ntuples, no one of which is linearly dependent on the others.
The dimensional anomalies, expressed in the theorems above, lead CuninghameGreen to introduce the idea of strong linear independence.
38
Definition. Vectors vi, ... , E W" are called strongly linearly independent (SLI) if and only if there exists a finite element of W', i.e. a vector v E V" such that v has a unique representation
V = I (ai * v') , ai E WIf V C Wnxk is the matrix whose columns are v, ... , Vk C W", then the strong linear independence of vi, ... , vk E W" is equivalent to the following condition: There exists a vector v C V' such that the equation V * x = v is uniquely soluble. Another equivalent condition is given by the theorem below, whose proof takes advantage of the following lemma:
Lemma 3.6.4. Suppose v1, ... , Vk E V. The vector v = V (nvi) * v E
is bounded from above by the ndimensional vector, whose entries are all equal to lv, i.e.: vj 5 IV V J = 1,... , n .
Proof. Let us consider the algebraic expression given by O_1 [(*"_1v) * 1v Opening with respect to the index j yields the following inequality, which is valid in V:
=4(c~1v1= *v0 07= I(~~ r * >(V)1 Vj = 1..I
< =iv VJ = I ... n
=lv V J = 1,...,n .
Theorem 3.6.5. The vectors vi, ... , vk E V" are SLI if and only if ok is
where the symbol i denotes the vector (Ov) * v for all v E V.
39
Proof. "SLI => Eq. (35)": If given vectors vi, ... ,vk C V' are SLI, then v = Z =1c E V" is the unique representation of this vector in terms of vi, ... , vk. Hence, k91 $ c19i for all p = 1,...,k. The inequalities O i for all
p = 1,. .., k imply that Eq. (35) holds.
"Eq. (35) > SLI ": Given vectors vi, ... , Vk E V', we assume that Eq. (35) holds, i.e. 71i" < =1i0 for all p = 1,..., k. We set ai = (1vi>1= (Ov')1. We will show that the vector v = O (a, * vi) C Vn does not have another representation. Hence, if v = K (ci * v'), it suffices to show that c= ai for all i = 1,..., k.
The following argumentation shows that ci < ai for all i = 1,..., k. The statement v = Kf_ (c, * v') is equivalent to v = i (ci * V') for all j 1, . . . , n. For fixed, but arbitrary j E {1, . . ., n}, consider the algebraic expression K4 1 (c * oj). An application of the principle of opening with respect to the index i yields:
ci *o V 5 1(i V Vi 1..k .
By assumption, O_ (ci * v) = vj for all j = 1,... , n, and by the previous lemma the entries vj of v are bounded from above by Iv. Therefore,
c,*v < lv Vi=1,...,k; Vj=1,...,n .
A right multiplication by (v;) yields:
ci < (v') Vzi = 1,..., k; Vj= 1,...,n .
We are able to finish the proof of ci < ai for all i = 1,... , k by closing with respect
40
to the index j and by using Eq. (32):
C1. < (V [K (v')K = (Ovyl =i
Assume that there exists a p E {1,..., n} such that cp < ap. Then c, * v < ap * v for all j = 1,... ,n. By the principle of opening, ai* v< O_=a * v ) = vi V i =1,...,k; V j=1,...,n . Thus, cp *v < vJ for all j = 1,...,n, which implies that c * vP < v. Combining the inequalities Q 1(ci * v') < v and c, * vP < v, where j = 1, ... ,n, leads to a contradiction to the assumption v = K{1(ci * vi): Q (ci * v2) = (Ci * V') 0 (Cp * v ) < v .
Example. Let the following elements of R3 be given: Vi = ,0 V2 = ,' V3 = 2 .
Note that
 the corresponding normed vectors i' , i = 1, 2, 3 are given by i 1=( 3 ,2= 0 ,#3= 5 .
0 3 2)
" The vectors vi, v2, V3 are not SLI since 9l V #r2 = ( = fi .
(0 i=1
41
The following theorem shows that the concept of strong linear independence is more suitable for defining notions of rank and dimension in minimax algebra than the concept of linear independence.
Theorem 3.6.6. For any linear blog W there are k vectors vi... ,Vk E W" which are SLI if and only if k G {1,...,n}.
Definition. We define the SLIrank of a finite matrix A E Vmxn, denoted by ranksLI(A), as the maximal number of row vectors which are SLI. Theorem 3.6.7. The SLIrank of A E Vmxn equals the maximal number of column vectors of A which are SLI.
Note that by Theorem 3.6.6 the SLIrank of a finite matrix A E Vmxn is less than or equal to min {m, n}.
In linear algebra, the concept of matrix rank naturally induces a concept of rank of rectangular templates. This notion has found many applications in the area of template decomposition due to the fact that it represents the minimal number of those products of row and column templates which add up to the original template. These considerations lead us to introduce a concept of matrix rank within minimax algebra which has similar consequences with respect to the purpose of template decomposition. Definition. We define the rank of a finite matrix A E Wmxn as the minimal number of vector pairs into which it can be decomposed, i.e. the minimal number r for which there exist column vectors ul, ... ,ur E Wmx and row vectors w', ... wr E Wl"", such that A has a representation of the following form: A =41 * VI.
42
Theorem 3.6.8. Let W be a linear blog. The rank of a finite matrix A E W' is bounded from below by the SLIrank of A and bounded from above by the maximal number t of linearly independent row vectors of A as well as by the maximal number of linearly independent column vectors of A. Proof. "ranksLI(A) < rank(A) ": Recall that the rank of A is defined as the maximum number k of SLI row vectors. Let v, ... ,vk E W"' be k SLI row vectors of A. By Theorem 3.6.5
Suppose that rank(A) = r. Then there exist row vectors d.,... d' E W1"' such that each row vector of A can be expressed as a linear combination of d,....., d'. In particular, each row vector v , where i = 1,... , k, has such a representation. Hence, there exist scalars c, C W, such that the normalized row vectors v' have the following representation:
1 = c *d Vi=1,...,k.
Let cj = _ (c') and let v = _ ir. Using the axioms of associativity, commutativity, and distributivity, we show that v = (c, * d): V = Ok'01=1 c0_1 (c' * d)]
= j=1 [Or1 (c' * d')]
i(, (c)) * d
= c ;( * d').
43
Given any I E {1,. . ., r}, there exists and index i {1,C . . , k} such that cl = c' due to the linearity of W. This fact yields cl * dl < v' in light of the equality Or c * dl) = v'. Closing this inequality with respect to the index 1 yields r (c; * d') O{_1v'. Thus, we are able to conclude the proof of the inequality rank(A) > ranksLI(A) as follows:
V=O (c, *d') 1= iv  S=,v~
. 4 v
4 ii ..,i,}> k
r>k.
"rank(A) < t ": It suffices to show that there exist column vectors ul, ... ,ut E W 'Xl and row vectors w', ... , wt E Wl"', such that A has a representation of the following form:
A = 1 ( M v).
Let a(i) E Wx"' denote the ith row vector of A for all i = 1,...,m. WLOG the first t row vectors of A are linearly independent, so each row vector a(i) can be expressed as follows in terms of a(1),...,a(t): a(i) = 01'=(c' * a(i)) Vi = 1,..., M ,
where cj E R are some appropriate scalars for all i = 1,..., m and for all t = 1,..., 1. If cl E R"'"x denotes the column vector
Cl = (') V ,.,t
c,
44
we have A = (1 c' * a(l)
which concludes the proof of the theorem.
CHAPTER 4
RANKBASED DECOMPOSITIONS OF MORPHOLOGICAL TEMPLATES
4.1 Introduction
By way of a canonical bijection, matrices correspond to rectangular templates. Hence, the notion of matrix rank induces a notion of rank of rectangular templates both in the linear domain and in the nonlinear domain. Moreover, algorithms for the decomposition of matrices translate into algorithms for the decomposition of templates. A rankbased decomposition of a rectangular template is considered to be a weak decomposition of a rectangular template according to its rank, i.e. a decomposition into a minimal number of sums of outer products of row and column templates. Rankbased decompositions have the following benefits:
0 The computational complexity of a convolution of an image with a rectangular
template reduces from a problem which is quadratic in complexity to a linear problem. 0 The rankbased decomposition of a template can be used as a first step for the
decomposition into a form which matches the neighborhood configuration of a
particular parallel architecture.
The strong decomposition of a rank I template is an easy task both in the linear and in the nonlinear domain [38]. O'Leary [44] showed that any linear template of rank 1 can be factored exactly into a product of 3 x 3 linear templates. Templates of higher rank are usually not as efficiently decomposable. However, the LU factorization yields a proven method for determining a rankbased decomposition of a linear template of arbitrary, unknown rank [28, 48].
45
46
Thus far, rankbased decompositions of morphological templates are restricted to decompositions of rank 1 templates. In this chapter, we provide a polynomial time algorithm for the rankbased decomposition of morphological templates of rank 2. Furthermore, we present a heuristic algorithm for the rankbased decomposition of morphological templates of arbitrary rank.
4.2 Correspondence between Matrices and Rectangular Templates
Note that there is a natural bijection 0 from the space of all m x n matrices over R,. into the space of all rectangular m x n templates in (RX,0)X.
Let y = (X, y) E X be arbitrary and x = (x1, x2) E X be such that
X1=xmn 2M 1 , X2=ymin , n I
The image of a matrix A E R""" under 0 is defined to be the template t E (R?,)X which satisfies
ty(X1 + i  1, X2+ )= aij V I =,.,m, V j ,.n,
ty(yi, Y2) = 00 V Y1, Y2 [X1, X1 + M  ]X [X2, X2 + n  1].
The natural correspondence between rectangular templates in t E (Rio)X and matrices over R, allows us to use a minimax algebra approach in order to study the weak decomposability of rectangular templates into separable templates.
Example. Let A E R3x3 be the matrix and u, v the vectors given below.
2 510
A= 4 7 12 , = 0 v = (7,4,12).
3 0 5 7
47
The function 0 maps A to the square template r C (R,)X in Figure 1, and it maps the column vector u to the column template s E (Rxo)X and the row vector v to the row template t E (RX,0)X in Figure 6.
Definition. If t = O(A) for some realvalued matrix A, then we define the rank of the template t E (Ro)X to be the rank of the matrix A.
Our interest in the rank of a morphological template is motivated by the problem of morphological template decomposition since the rank of a morphological template t E (RX.)X represents the minimal number of separable templates whose maximum is t or, equivalently, the minimal number r of column templates r' E (Rx.) X and row templates s C (Rx..) such that t = V (rI M si).
i=1
Lemma 4.2.1. If A E Rx"' and u Mv < A, where u C R"x1 and v E RX"', then u < A (aj  v) for all i = 1,...,M.
j=1
Proof. Each entry aij of A is greater than or equal to ui + vj for all I = 1,..., k,
n
implying that ui A A (ai2v ) for all i = 1,... M.
j=1
k
Theorem 4.2.2. If a matrix A G R mx has a representation A = V (u' M v') in terms l=1
of column vectors uc E R"' and row vectors vi E Rxn, where I = 1,..., k, then A can be expressed in the following form:
k
A = V (w'Mv'
l=1
where w1 E R' is given by
n
w! A (aij  v) (i=1,...,m).
j=1
48
k
Proof. Let W = \/ (w' M v1). The proof is accomplished in two parts.
1=1
1. "W < A":
k
=1
k
V W(+ v. 1=1
I .j
I
= a ,j
z=1,...,m, j=1,...,n.
k
2. "W > A": Since A C R"" can be represented as V (u' M v1), we observe that 1=1
by the previous lemma u (aij  v1) for all i, j, and 1. The latter inequality has
the following consequences, which conclude the proof of the theorem:
k
a,3 = V (u! +V
1=1
< k m ~, k
V Q R (ajs v + v = V wij + .
Remark. Theorem 4.2.2 implies that, for any matrix A E R"" of separable rank k, it suffices to know the row vectors vi E RX"<, I = 1,...,k, which permit a weak decomposition of A into k separable matrices in order to determine a representation of
k
Wi =1M ij
49
A in the form
k
A =V (W' MV') l=1
w1 E Rm' V l= 1,...,k. Like most of the theorems established in this dissertation, Theorem 4.2.2 has an obvious dual in terms of column vectors which we choose to omit.
4.3 Invariance Properties of the Rank of Realvalued Matrices
In this section, we derive some invariance properties of the matrix rank which translate directly into results about the rank of rectangular templates. These theorems greatly simplify the proofs of the decomposition results which we will present in the following sections. Definition. Let A E R' " and p be a permutation of {1,..., n }. The associated column permuted matrix pc(A) of A with respect to p is defined as follows:
alp(,) alp(2) ... alp(n) a2p(l) a2p(2) ... a2p(n) pc(A)= amp(1) amp(2) ... amp(n) Similarly, if p is a permutation of {1,..., m}, then we define the associated row permuted matrix p,(A) of A with respect to p as follows: a,(1)l ap(1)2 ... ap(1)n Pr (A) = ap( :2)1 ap( 2)2 ... ap( .2)n ap(m)1 ap(m)2 ... /p(m)n
Theorem 4.3.1. The rank of a matrix A e R"' is invariant under column and row permutations.
50
Proof. Suppose the matrix A E R"'" has rank r. Hence, there exist column vectors u1 C R"' and row vectors vt E R"', such that A can be written in the form
r
A = (ul E VI),
1=1
Let p be a permutation of the set {1,...,n}. An application of the column permutation pc to the matrix A yields the following identities for all i = 1,... ,m and for all j = 1,...,n:
r r /r
(pc(A))ij = agpg) = + V(j) = V (U' M (= U1 M Pc(V)
1=1 1=1 =i
Hence, pc(A) =V ul M pc(vl) , which implies that the rank of pc(A) is at most r.
1=1
Now it suffices to show that pc(A) can not be weakly decomposed in the form
k ~
pc(A) = V ul M v1 , where ul E Rmixi, v' E R,)xn, and k < r. Suppose, we are given
l=1
k~
a representation pc(A) = V ul M vI . Applying the inverse of the column permutation 1=1
pc to both sides, we obtain the following equation:
k
A = \ 1 to pC 1 .
l=1
This equation concludes the proof of the invariance of the matrix rank under column permutations, because k must be at least as big as r = rank(A). The invariance of the matrix rank under row permutations can be proven in a similar fashion.
Recall from elementary linear algebra that the set of realvalued m x n matrices together with the operations of addition and scalar multiplication forms a vector space over R. The operation of addition in R"' < is performed pointwise and denoted as usual by the symbol +. This addition is not to be confused with the operation V, the addition
51
in the band space R"'. Likewise, we differentiate between the scalar multiplications of the vectorspace R"' " and the scalar multiplications of the band space R"x".
Theorem 4.3.2. The following operations of the vectorspace R"' preserve the rank of a matrix A C R"' in minimax algebra: a. additions of matrices of rank 1;
b. scalar multiplications.
Proof. Let r be the rank of the matrix A E R"". By definition of the matrix rank, the integer r is the minimal number such that A V (u' M v') for some column vectors 1=1
u1 E R" < and row vectors vI E R',n. Suppose, that we are given a separable matrix B C R"'x and a scalar c E R. The separability of B implies that B can be written as an outer product of a column vector u E Rmx, and row vectors v E Rx.. The validity of the theorem follows from the following facts, which are easily derived: A+B= \ u + Uu) M (v + V)]
c  A =V (c  ul M c  VI.
1=1
4.4 Separable Matrices and Templates
At this point, we are finally ready to tackle the problem of determining rankbased decompositions of matrices in minimax algebra. The reader should bear in mind the consequences for the corresponding rectangular templates.
For any matrix A E R"x", we use the notation a(i), where i = 1,..., m, to denote the ith row vector of A and we use the notation a[j], where j = 1,.. . , n, to denote the jth column vector of A.
52
Theorem 4.4.1. (Li and Ritter [38]) Let A E R"'xn be a separable matrix and { a(i) : i = 1,... , m} be the collection of row vectors of A. For each arbitrary row vector a(i0), where 1 < io m, there exist scalars Ai E R, i = 1,... ,m, such that
the following equations are satisfied:
Ai +a(io,) = a(i), V I=,...m.
In other words, given an arbitrary index 1 < io 5 m, each row vector a(i) is linearly
dependent on the ioth row vector a(i0) of A.
Corollary 4.4.2. If A e R'" is a separable matrix then each column vector of A is linearly dependent on the j,th column vector, where j' is a fixed, but arbitrary column
index 1 < jo < n.
Proof. The proof follows from the fact that the separability of A E R' is equivalent
to the separability of AT E R n", the transpose of A.
Clearly, Theorem 4.4.1 gives rise to the following straightforward algorithm which
tests if a given matrix over R is separable. In the separable case, the algorithm computes
a vector pair into which the given matrix can be decomposed.
Algorithm 1. Let A E R" n be given and let a(i) denote the ith row vector of A.
The algorithm proceeds as follows for all i = 2,... , rn: STEP 1:
Form ci E R, where ci = ail  a1,1. STEP 2(j  1), j = 2,...,n:
Subtract ai from aij.
53
STEP 2(j  1) + 1, J = 2,...,n:
Compare ci with aij  aij. If ci aij  aij, the matrix A is not separable and the algorithm stops. If step 2(n  1) + 1 has been successfully completed, then A
is separable and A is given by c M a(1).
Note that this algorithm only involves (m  1)n subtractions and (m  1)(n  1) comparisons. Hence, the number of operations adds up to 2(m  1)n  m + 1, which implies
that the algorithm has order O(2mn).
As mentioned earlier, the given image processing hardware often calls for the
decomposition of a given template into 3 x 3 templates. In the case of a separable square template, this problem reduces to the problem of decomposing a column template
into 3 x 1 templates as well as decomposing a row template into 1 x 3 templates.
Theorem 4.4.3. Let t be a square morphological template of rank 1, given by t = r M s, where r is a column template and s is a row template. The template t is decomposable into 3 x 3 templates if and only if r is decomposable into 3 x 1 templates and s is decomposable
into 1 x 3 templates.
Proof. "4=": Let r', . . ., rk be 3 x 1 templates such that r = r' M ... M rk, and s, ... , sk be 1 x 3 templates such that s = s' M ... M sk. The template t decomposes
strongly into the 3 x 3 templates r' M s', z = 1, . . , k, for the following reason: t = r M s
=(r' . M rk) to (S . sk)
= (r' M s' ) M ... M ( rk M Sk) .
54
">": Let t = r Ms be strongly decomposable into k 3 x 3 templates. We show by induction on the number k that r is decomposable into 3 x 1 templates and s is decomposable into 1 x 3 templates.
If k = 1 then t is a 3 x 3 template and therefore, r is a 3 x 1 template and s is a 1 x 3 template which leaves nothing to show. For the remainder of the proof, we assume that if a square template t = r M s, where r is a vertical and s is a horizontal strip template, decomposes into less than k 3 x 3 templates, then r is strongly decomposable into 3 x 1 templates and s is strongly decomposable into 1 x 3 templates.
Let a template t of rank 1 be given which can be represented in terms of 3 x 3 templates as, i = 1,..., k.
t = a' M ... M a*
Let us define the template a as a' M ... M aki, and the template b as a . Furthermore, we introduce the following convenient notations for any rectangular p x q template. Recall that for any matrix C E RPXq the symbols c(i) and c[j] denote the ith row vector of C, the jth column vector of C respectively. The indices i range from 1 to p and the indices j range from 1 to q. Let us introduce the following convenient notations for arbitrary rectangular templates c.
c(i) = (1(c)(i)) C( (41(c(i))) ,
c[j = (1(c)[j) Cu = (1(c[j]))
Obviously, t = a M b. Also note that the row template t(1) and the column template t[1] can be computed as follows:
t(1) = a(l) M b(l) , t[1] = a[1] M b[1] .
55
Since t is of rank 1, all row vectors are linearly dependent on q~1(t(1)).
Hence, there exist constants ci, i = 1,..., m such that t(i) = ci + t(1). In fact, ci = t(i),  t(1), for all i = 1,... , m, which implies that t has another representation besides t r M s, namely t = (t[1]  t(1)1) M t(1). We conclude that ri+s = t[1];+t(1) for all i = 1,... , m, in particular for i = 1. Thus, s and t(1) only differ by a constant c, and for similar reasons r and t[1] differ by a constant d.
s t(1) + c
r = t[1] + d.
We are now able to finish the proof by substituting a(1) M b(1) into t(1) and a[1] M b[1] into t[1].
S =W() + c =[a(l) M b(1)] + c =a(l) M [b(1) + c]
r t[1] + d (a[1] M b[1]) + d = a[1] M (b[1] + d)
Example. Let A be the realvalued 5 x 5 matrix given below.
3 7 4 2 1 11 15 12 10 9 A= 5 9 6 4 3 = uMv,
6 10 7 5 4 \12 16 13 11 10 J where
/0
8
u = 2 , = (3 7 4 2 1).
3
The template t =O(A) is not decomposable into two 3 x 3 templates since the template r = q(u) is not decomposable into two 3 x 1 templates.
56
4.5 Matrices and Templates having Rank 2
Lemma 4.5.1. If a matrix A E R"" has rank 2, then there exits a transform T consisting of only row permutations, column permutations and additions of separable matrices  as well as vectors u E R'" 1 and v E R Xfl such that T(A) can be written in the following form:
T(A) = Omxn V (u v) ,
where
U1 < U2 <   <. Um ,
and 0.,x , denotes the m x n zero matrix. Proof. The matrix A E R" " can be represented by means of some column vectors u', u2 E R"'X ' and some row vectors v', v2 E Rlxn: A= (u'Mv1) V (u2Mv2). If B = (ul) M (v') then we compute the matrix A + B E R" " as follows: A + B = [(u' M v') V (u2 M v2)] + B = [(u' M v1) + B] V [(u2 M v2) + B] = Omxn V [(u2 _ u') M (v2  ,] Let a denote the column vector u2  u' and let b denote the column vector v2  v'. The entries of this vectors can be ordered by a permutation p : {1,... , m}  {1,... , m}
57
such that
ap(1) < ap(2) < ... ap(m)
Applying the row permutation pT to A + B yields the following result: p,(A + B) =
Pr(0mxn V [a M b]) =
pr (Omxn) V pr (a M b) =
Omxn V [pr(a) M b].
Hence, if we introduce the notations u for pr(a), v for b and T for the addition of B followed by pr, we are able to conclude the proof of Lemma 1: T(A) =0 V (uMv),
U1 U2 < ... < Um .
Lemma 4.5.2. Let T one of the rank preserving transform of Theorems 4.3.1 and 4.3.2 and A E R'"*". The transform T maps row vectors of A to row vectors of T(A), and column vectors of A to column vectors of T(A). Proof. Pick two arbitrary entries aik and ail on the ith row of A and observe where they are mapped under row permutations, column permutations, additions of row vectors, additions of column vectors, and scalar multiplications. Theorem 4.5.3. A matrix A E Rmx' has rank 2 if and only there are two row vectors of A on which all other row vectors depend (linearly).
58
Proof. Since A C R"' has separable rank 2 there are vectors u', u2 E R"'1 and vectors v', v2 E Rx"' such that A can be expressed as follows: A = (u' M v') V (u2 M v2).
By Lemma 4.5.1, there is a transform T and vectors u E R"' and v E R " such that T(A) = Omxn V (u M v), U U2 <    < Um .
Let B denote T(A). Note that for every index i {1,... ,m} the ith row vector b(i) is given by b(i) = 0x V (ui + v). The ascending order of the scalars ui, z = 1,... ,m, induces an ascending order of the vectors b(i), i 1,...,m. Using the fact that u1 < u2 :5 ... urn, we can also show that b(i) is linearly dependent on b(1) and b(m):
b(i) = O1,n V (ui + v) = Oixn V [(ul + v) V (ui + v)] = Ox V [(ul + v) V (u2  UM + UM + v)] = [Olx, V (ui + v)] V [(ui  un) + (Um + v)] = b(1) V [(ui  urn) + b(m)] . These relationships between the vectors b(i) and the vectors b(1) and b(n) imply that B can be written as follows:
B = [0.x M b(1)] V [(u  urn) M b(m)] .
59
An application of the inverse transform T1 yields the following: A = T1(B) =
[T1 (O.xi) M T1(b(1))] V [T'(u  um) M T1(b(m))]
Recall that, by Lemma 4.5.2, the vectors T1(b(1)) and T1(b(m)) are row vectors of A. Hence, there are indices k, 1 C {1,...,rm} and vectors c', c2 E R"' ', such that A = [cl M a(k)] V [c2 a(l)]. In other words, every row vector a(i), where i E {1,... , m}, is linearly dependent on the same pair of row vectors a(k) and a(l). Therefore, we have completed the proof of the sufficiency part of the theorem.
Obviously, the other direction of the theorem holds.
Remark. By Theorem 4.5.3, a matrix A E R"'f has rank 2 if and only if there exist two row vectors of A  a(o), a(p) where o, p E { 1, ... , m}  which allow for a weak decomposition of A. In this case, an application of Theorem 4.2.2 yields the following representation of A:
A = [u M a(o)] V [v M a(p)],
where
U, v E R"x,
Ui = A (aaoj ) (i= 1,...,m),
In
vi = A(aij  aj) (i= 1,...,m) .
=1
Hence, in order to test an arbitrary matrix A E R"x" for weak decomposability into two vector pairs, it is enough to compare A with (b[i] M a(i)) V (b[j] M a(j)) for all
60
indices i, j E {1,..., m}, where the matrix B E R"m"' is computed as follows:
n
bi = (a,8  aj,) 8=1
vi i= 1, ...,mI VJ = 1, ..., .
Algorithm 2. Assume a matrix A E R"x needs to be decomposed into a weak M product of 2 vector pairs if such a decomposition is possible. Considering the preceding remarks, we are able to give a polynomial time algorithm for solving this problem. For each step, we include the number of operations involved in square brackets. STEP 1:
Compute B E R"'.
[m2n subtractions; m2(n  1) comparisons] STEP 2:
Form C' = b[i]M a(i) for all i = I,., m.
[m(mn) additions] STEP 3:
Form C V C' for all i, j =1,.. ., m and compare the result with A. If C" V CP = A
for certain o, p = 1,... , m then the algorithm stops yielding the following result: A = (b[o] M a(o)) V (b[p] M a(p)) .
If C' V Ci < A for all i, j = 1,..., Im, then A does not have a weak decomposition
into two vector pairs.
[At most 2 ) (mn) comparisons]
This algorithm involves at most a total number of m2(3n  1 + (m  1)n) operations
which amounts to order O(m3n).
61
Example. Let us apply Algorithm 2 to the following matrix A E R45
9 5 5 6 5 A= 7 4 4 3 4
12 7 7 9 7 11 3 5 8 3 We compute matrices B E 11414 and C' E 11415 for all I=1,.,4
0
B= 3
2 2
C' =b[1] M a(1) = C2 =b[2] M a(2) = C3 =b[3] M a(3)
C4 = b[4] M a(4) =
1
0
3
1
3
6
0
4
9 C6
11
7
8
7 10
6
9 62
7
9 C6
12 11
2
5
1 '
0 /
5
2
7
3
5
4
7
3
4
1
7
3
1
2
4
3
5
41 7'
3
5
4 7'
3
4
7 ' 3)
1
2
3
Comparing the matrices C V Ci with A for all i, j = 1,..., m reveals that A = C2V C4. Thus,
A = (b[2] M a(2)) V (b[4] M a(4)) .
Example. Let A E R9x9 be the following matrix, which constitutes the maximum of a
matrix in pyramid form
18
14 10
6
A= 2
6
10
14
\18
and a
14
9
4
1
2
1
4
9
14
matrix
10
4
1
4
6
1
1
4
10
62
in paraboloid form.
6 2 6 10
1 2 1 4
4 6 4 1 7 10 7 4 10 14 10 6 4 6 7 4 4 6 4 1
1 2 1 4
6 2 6 10
Since matrices in paraboloid and in pyramid form are separable, algorithm 2 should yield a weak decomposition of A in the form [u M a(o)] V [v M a(p)] for some vectors u, v E R9 and some indices o, p E {1, .. ,9}. Indeed, if B E R9X9 denotes the matrix computed by algorithm 2,
A = [b[1] M a(1)] V [b[3] M a(3)],
where
b[1] =
a(1) = (18
14 
0
4 8
12 16 ,b[3]=
12
8
4
\0 10 6 2
11
5
0 3 10 3 0 5
\11/
6 10
14 18) ,
a(3)=(10 4 1 4 6 4 1 4
9).
This weak decomposition of A can be used to further decompose the square templates O(b[1] M a(1)) and 0(b[3] M a(3)) into 3 x 3 templates. By Theorem 4.5.3, the templates k(b[1] M a(1)) and 0(b[3] M a(3)) are decomposable into 3 x 3 templates if and only if the
14
9
4
1
2
1
4
9
14
18
14
10
6
2
6
10
14
18/
63
column templates 4(b[1]) and 4(b[3]) are decomposable into 3 x 1 templates and the row templates q(a(1)) and 4(a(3)) are decomposable into 1 x 3 templates. It is fairly easy to choose 3 x I templates ri E (Ri0)X and 1 x 3 templates si E (R3,) X,i = 1, ..., 4, such that 4(b[1]) = ((r1 M r2) M r3) M r4 and 4(a(1)) = ((s1 M s2) M s3) Ms4. For more complicated examples, we recommend using one of the integer programming approaches suggested in [23, 58, 57].
r1 = r2= r3= r4= y y y y
0
4
0
s = 18 14 18 s 2= s 3= s= 0 4 0
Figure 7. Templates of size 3 x 1 and size 1 x 3
providing a decomposition of the templateO(b[1] M a(1)). Hence, we obtain a representation of O(b[1] M a(1)) in the form below.
q(b[1] M a(1)) = [((r1 M r2) M r3) M r4] M [((sI M s2) M s3) M s4]
By rearranging the templates ri and s' for i = 1,... , 4, we can achieve a decomposition of O(b[1] M a(1)) into four 3 x 3 templates t' = r' M si.
y
Figure 8. 3 x
18 14 18
14 10 14 t2= t3 t=
18 14 18
3 templates providing a decomposition of the
0 4 0 4 8 4 0 4 0 template 0(b[1] M a(1)).
64
In a similar fashion, we are able to decompose the template 4(b[3] M a(3)) into four
3 x 3 templates.
4.6 Matrices and Templates of Arbitrary Rank Unfortunately, rankbased decompositions of matrices having rank k > 3 can not be obtained by employing a similar strategy as in Algorithm 2, since a theorem similar to Theorem 4.5.3 does not hold for this class of matrices. This fact is expressed in the following theorem.
Theorem 4.6.1. For every natural number k > 3 there are realvalued matrices which are weakly M decomposable into a product of k vector pairs, but not all of whose row vectors are linearly dependent on a single ktuple of their row vectors.
Proof. Let n = ( for any given k C N. In order to prove this theorem we show
that a matrix A E Rnxk satisfies the statement of the theorem if it has the following properties: Each row vector a(i) of A has two locations i,, i2 E 1,. . . , k} having value
1. All other entries of a(i) are 0. If i j, then a(i) / a(j).
Note that a(i) = v V vi: for alli = 1, ...,n, where v E R I1 X is defined as follows for all I = 1,.. k:
0 else .
k
Equivalently, the matrix A has representation A = V (ul M v'), where the column l=1
vectors ul E Rnx 1 (1 = 1,.. ., k) are determined by
_i I if I = i1, Z2
1 0 else .
65
Therefore, A satisfies the first requirement of the theorem.
Let (a',. .., ak) denote an arbitrary ktuple of row vectors of A. (Note that the index i in a' does not necessarily correspond to the row index.) Assuming this ktuple allows for a weak decomposition of A into k vector pairs, there exist column vectors
k
b E R nx (1 = 1,..., k) such that A can be written as V (bl M a). In other words,
k
a(i) =V (b + a' ) V I =1.. n . (41)
1=1
Now let us pick an index i E {1,. . ., n} such that a(i) $ a' for all I = 1,... k. For notational convenience, let us use the abbreviations a for a(i) and bl instead of b( for all i = 1,...,n and for all I = 1,...,k.
Since a differs from every vector a among al,..., a , there has to be an index 1' for each I = 1,..., k such that al, = 0 and a', = 1. Since bl + a' < a, the scalar bl cannot exceed 1. These considerations lead to the following sequence of inequalities contradicting the assumption presented in Eq. (41): k k
V (b' + a') V V[(1) + xk = Oxk < a. 1=1 1=1
Hence, the row vector a of A is linearly independent from the ktuple (al,... , ak).
Corollary 4.6.2. Let k > 3. There exist integers m and n, and matrices A, B E R'"** such that rank(B) = k and such that for any metric d on Rm x"
d(A,Vuita(i) >d(A,B)
VIC{1,...,n} I=k;
VUiER'
66
Proof. The previous theorem implies that for every k > 3, there exists a realvalued matrix A satisfying
d A, \V (u" M a(l)) > 0 = d(A, A)
iEI /
VIC j 1,...,n} : III = k;
V ui E Rm
We will show in chapter 6 that the general problem of determining the weak decomposition of matrices is NPcomplete. However, the NPcompleteness of a certain problem does not preclude the existence of an efficient algorithm for solving arbitrary instances of this problem. For matrices having relatively small rank compared to their size, we suggest the following heuristic algorithm. This heuristic is based on the following observations:
Remark. Suppose that A C Rmx", whose rank is an unknown integer r. Hence, A can be represented as the maximum of matrices A' E Rmx", where I = 1,...,r and A' = x1 M y' for some realvalued column vectors xi of length m and some realvalued row vectors yi of length n. According to Theorem 4.4.1, the separability of the matrix A' induces the equality of the differences a  a(2 for all i = 1, . .. , m and for arbitrary, but fixed 1 E {1,...,r} and j1,J2 E {1,... ,n}. On the other hand, if ai3,  aij2 for i E I C {1,...,m}, it seems reasonable to assume that ai2, = a = x! + y1 and aij2 =ak= x2 + ya( for most i E I provided that "r = rank(A) is relatively small and that I{a3i : i = 1,..., m; j 1,... ,n} is relatively large compared to r". This assumption and the following two lemmas play an important role in Algorithm 3 which intends to determine r = rank(A) as well as column vectors ul E R"' i and row vectors
67
V1 E R1x" such that
r
A = (ul to V1).
l=1
Lemma 4.6.3. Let A E R"" and x,u C R"x 1 and y,v E R' n. Suppose that there exists an index j' C {1,... , n} such that xi + yj, = aij, = ui + v3' for all i C I C {1, ... , m}. If for all j E J C {1, . . . , n} there exists an index ij E I such that ui, + vj = ai,j xi, + yj, then
ui + vj > xi + yj V i E I, V j E J.
Proof.
ui + v3 = (ui  u.,) + (ui, + V3)
= [(ui + vj')  (u,1 + Vji)] + (u2, + Vj)
= [aiji  aij] + (ui, + vj)
S[(Rxi + y)  (xi, + y)] + (xi, + yj)
= [xi  xi,] + (xi, + yj)
= xi + yj
V i C I, V J E J.
Lemma 4.6.4. Suppose u, x E R",x1 v,y E R'x2, and ji,j2 E {1, . ., m}. Furthermore, let a subset I of{1,... , m} and a subset J of {1,... , n } be given such that v3 and xj are finite for all j E J. If the following conditions are satisfied then ui + v3 xi + yj
68
for all i = 1,...,m and for all j C J.
ui+vj, =xi+y1j = ai,1 ViE I, ui+vj2 = xi+yj2 = aij2 ViE I, ui +v, = aij, or u, +vj,2 = aij2Vi' E {1,...,m} \ I, ui+vj >xi+yj ViE I; VjE eJ.
Proof. It suffices to show that for arbitrary i' E {1,... , m} \ I and for arbitrary j E J the relation ui' + v) > xi, + yj holds. Suppose j' denotes Ji if ui' + vj, = aij, and J2 otherwise. We obtain the following estimate for the difference ui  ui,, where i is an arbitrary element of I: U  Ui = (uj' + vj')  (u, + vji) = agl  aip
(xi + y+ )  (Xi + y) = xi,  xi
Using this estimate, the inequality ui, + vj xi, + yJ now follows from the inequality ui + V3 Xi + yj.
k
Remark. Assume that A V (x' f y'), where x' and yi are unknown realvalued
1=1
vectors. Algorithm 3 constructs vectors ul and v1 for I = 1,... , r such that ul M v1 < A.
k k
The maximum \/ (u' M v') reaches A, i.e. V (ul M v1) = A, if there exist sets P1 and
1=1 1=1
69
J1 for all I = 1,..., k such that I x J C{(i, j) i + yj = aij} I x Ji {(i,j) : i + j < Ui + Vj}
k
J (Ix Ji) ={1,..., m} x {I,..., n}. 1=1
Lemmas 4.6.3 and 4.6.4 already provide some insight into these matters. After presenting Algorithm 3 and applying this algorithm to a specific example, we establish a theorem which should lead to an even better formal understanding of Algorithm 3 by combining Lemmas 4.6.3 and 4.6.4.
Algorithm 3. Let A E Rm"'. This algorithm determines a parameter k > rank(A) and progressively constructs column vectors ul C Rx1 and row vectors vi E Rlxfl such that
k
A= V (ulMV)
1=1
Furthermore, this algorithm generates a matrix M e Rmx" which has the following property:
m = > u, + v = aij Vi = 1,..., m; Vj = 1,..., n; V l= 1,..., k
k
Informally speaking, each entry mij represents a marking indicating where \ (u + '4) is adopted.
Initialize the matrix M E R'xn to be the zero matrix Omxn and let 1 = 0. As long as the matrix M has a zero entry, we increment the parameter I by 1, and perform the following steps. Otherwise k = I and the algorithm stops.
70
STEP 1: (Find the most frequent tuple difference)
Determine a strictly increasing sequence of maximal length (i, i. ., iP) such that
there exist two different column indices J1 and J2 having the following property: asilj  aiJ2 = ai2jl  ai2j2 = . . =  aij2 and
mitil = 0 or Mitj2 = 0 Vt= l, ,p For any t = 1,...,p satisfying m,jl = misj2 = 0, we set mi,j, = rnij2 = 1. If either one of the entries mi,j, or Minj2 differs from 0 for all t = 1,... ,p, then set
mij, = rnij2 = 1 for all t = 1,...,p.
STEP 2 (Extend the markings of found tuples horizontally):
For all i {ii,. . .,i,} and all 1,... n, we set mij = I if
a. mi, = 1.
b. aij  aij, =A (azj  azj,).
z=1
STEP 3 (Create preliminary row vectors v1 E R"_l):
Initialize the vector v1 E R' by setting v = oo for all j 1,..., n. Let s = min {i E {ii,. .. , i,} : mi, = I}. For all j = 1,..., n such that there exists an index i E {i, ij} satisfying mij = 1, change the value of v to aij + aj,  ai,, where i E {i1,...,i,}. Note that the value of v1 does not depend on the choice of
an element i of the set {iM,...,i,} f {i 1,...,m :nm3 = l}.
STEP 4 (Extend horizontal markings vertically and create column vectors ul E Rmxl):
Define the column vector ul E R"Xl as follows:
= (aij V)
j=1
71
Note that u adopts a finite value for all i = 1,... , m, since A is a finite matrix and ij is finite. Increase the values of all zero entries of M to I whenever aij  V = u,. STEP 5 (Final horizontal extension of markings and creation of final vectors vi E Rmxl):
1
For allj = 1, ...,n such that V is finite, we equate v3 with v . Otherwise, we define
mas (aij  u!). Any remaining zero entry mij is set to I if ai, =Uj + V .
Example. Let us first illustrate how Algorithm 3 works with a simple example. Suppose
A E R3x4 is the following matrix:
A= 5
6
6
3
6
3
3
7
4
8
A matrix M E R3x4 is initialized as 03X4. Then we execute STEP I through STEP 5 for the parameter I = 1.
Since ai3  ai4 = 1 for all i = 1, 2, 3, we J2 = 4. The matrix M E R3x4 now has form 0 0 1 M= 0 0 1 (0 0 1
set (ii,... ip) = (1, 2,3), ji = 3, and
3 3
Considering A (ai  az3) and A (azi  az3), we recognize that
z=1 z=1
3
A (azj  az3) = (8  3) A (5  3) A
A (az2  a3) (6  3) A (3  3) A Z=1
Therefore
0
M= 0
(1
0
0
1
1
1
1
(6  7) = 1 = a31  a33 (6  7) = 1 = a32  a33
1)
72
Initially, all entries of the vector v1 are set to oo. Since mi3 = 1, the parameter s = 1. When executing STEP 3, these entries adopt the following values: V1 a31 + a13  a33 = 6 + 3  7 = 2. v a32 + a13  a33 = 6 + 3  7 = 2. = ai3 + a13  ai3 = a13 = 3 Vi = 1,2,3. V= ai4 + a13  ai3 = (ai4  ai3) + a13 = 1 + 3 = 4 Vi = 1,2, 3.
3
STEP 4 computes a vector u' E R3x1 as follows:
4
u1 = (a1  v =)(8  2) A (6  2) A (3  3) A (4  4) = 0 ,
3=1
= A (a2  o = (5  3) A (3  2) A (3  3) A (4  4) = 0 ,
=1
S= A(a3  1) = (6  2) A (6  2) A (7  3) A (8  4) = 4
3=1
Since the vector v is finite, we have v1 = v in STEP 5.
Now the parameter I changes to 2, since there are still some zero entries of M left. These zeros of M are eliminated in STEP 1, since a11  a12 = a21  a22 = 2 and M11 = M12 = n2l = 0. Thus,
2 2 1 1
M= (2 2 1 1).
1 1 1 1
Note that STEP 2 through STEP 5 effect no change in the matrix M, since M > 03X4STEP 3 and STEP 4 create the following vectors V2 E R1i4 and u2 E R3x1: V2 = (8 6 oc oo),
2 2
73
STEP 5 produces a row vector v2 E R1x4, where v? = v? for i = 1,2. The entries V2 and vi are computed as follows:
3
i=1
3
V3= A (ai4  Ut?) = (4  0) A (4  (3)) A (7  (2)) = 43
i=1
Since all entries of M are strictly greater than 0, Algorithm 3 finishes yielding the following result:
A= (u'Mv')V(u2Mv2)
2 3 4)] V (3 (8 6 3 4)
2
5 3 3 4
=8 6 3 4).
6 6 7 8
Theorem 4.6.5. Suppose A > xIy, where A C R" , x E R"1, and y E R1"'. Let
u, iv, v denote the vectors uV, v1 constructed from the matrix A in Algorithm 3 for an arbitrary, but fixed parameter 1. Let j1, J2 be as in Algorithm 3. If I and J denote the sets
I{i' ,..m: u + jj=xi + yjl = ai,and uj+ 2 =Xi + yj2 =aij21
J= = 1,...,n : aij  aij, =A (azi  azj,) for some i E I
z=1
then
aij ui!j +Vj = u + > xi +yj Vi =,...,m; V j E J .
In particular, if xi+yj = ai3 for some i E {1, . . ,m} and some j E J, then ui+v, = ai3.
= 0 M(2
4
74
Proof. Note that I C {ii,..., iI,} since ai23  aij2 =  j2 for all i E I. Let us first show that aij u2 + j xi + yj for all i E I and for all j E J. By Lemma 4.6.3, it suffices to show that for all j E J there exists an index u,, + i3, = ai,j > xi, + yj.
By the definition of J, each index J E J entails an index ij E I such that mij is set to I in STEP 2 of Algorithm 3. If s denotes min{i E {i1,...,i,} : mn,3 = 11, then j = aij + a,,  ai,j, for all j E J. For an arbitrary, but fixed element j E J the inequalities i5, <; ai,j, + aj,  ai3j, hold for all j' = 1,.. . , n. Therefore, we are able to deduce the following equalities:
n
ui, +i, A (aij,  fyj') +
j'=1
= (aij  i3) + i'
= ai,j
V] E J.
Since A > x M y, the conditions of Lemma 4.6.3 are satisfied, yielding that ui + lj xi + yj for all i E I and for all j E J.
Provided that ui, + ij, = aij, or ui, + vj2 = aij2 for all i' E {1,... , n} \ I can be shown, an application of Lemma 4.6.4 concludes the proof of the theorem, since vj = vj for all j E J and A > u M 4v by the definition of u. Let us fix an arbitrary row index
n
Z' V I. By the definition of ui, as A (aij  ij), there exists a column index j' such j=1
that ui, = aij,  vj'. Since ui, is finite, ij is finite and equals vj'. Therefore, j' either equals J2 or is an element of the set J' which we define as the set of all column indices j such that ai3  aij, = A (as,  azj,) for some i E {ii,. . . ip}. Note that J' includes z=1
ji. The case I' = J2 leaves nothing to show. Assuming that j' C J', we finish the proof of the theorem as follows. We use the facts that A > u M ' and that aij, = ui + jl
75
for all i E {ii,...,ip}.
m
aiij,  a'3, / (azj,  azj,) = aij,  aij, for some i E {i1,.. . , p} z=1
> up + ij  ai, > (ui + pji)  (ui + ,) for some i E {i,. ..p =~Ui' + j a,3
Theorem 4.6.6. Algorithm 3 is a polynomial time algorithm. Proof. Let us analyze STEP 1 through STEP 5 for A C Rm" and for fixed 1 E N.
Since there are n (ni) different tuples of column indices ji and J2 m (2 22
subtractions are needed in order to compute the differences aij,  aij2 for all i = 1, . . ., m and for all j1,j2 = 1,...,n. For fixed i1 and J2, the comparison of the m real numbers aij,  aij,2 takes m 1) operations. Hence, these comparison operations are of computational complexity O(  n2M2). All other operations of STEP 1 and the operations of STEP 2 through STEP 5 have smaller computational complexity. Since the number of zero entries of M C R"x' decreases when performing STEP I through STEP 5, Algorithm 3 is a polynomial time algorithm. Theorem 4.6.7. Algorithm 3 produces the desired results. Proof. Suppose that Algorithm 3 is applied to a matrix A E R"x, yielding as an output an integer k, column vectors u E R"' and row vectors vi E Rlx" for all = 1, ... k, as well as a matrix M E Rn .
Since the final matrix M belongs to the set {1,... , k}mxn, it suffices to show that the following properties hold for all i, j, and 1:
I. u, + 1 < a,3.
76
2. mij =l u + v > aij.
n
The first property immediately follows from the fact that u equals A (aiz  vi), since z=1
n
A (aiz  vi) !; aij  V. z=1
Suppose that mij = I after completion of STEP 3. This assumption implies that i E {i1,.., iP} and v= = ai + a831  aijl. In this case, the following inequalities show that property 2 holds:
n
v A (aiz V +V
z=1
n
A (az  (aiz + aj,  ai3)) + V) z=1
= (aij,  a8,3) + (aij + aj,  aij)
= ai, .
Clearly u + = ai,, if the entry mij of the matrix M E R" " has been set to 1 in STEP 4 or in STEP 5.
Theorem 4.6.8. Given a matrix A E R'n of rank 1, Algorithm 3 determines a column vector ul E R"1l and a row vector v' e R1^, such that
A = u Mv1 .
Proof. By Li's Theorem, each row vector a(i) of A, where i = 2,..., n, is linearly dependent on the first row vector a(1). Hence, there exist Ai such that a(i) = Ai + a(1), which implies that the following equalities hold: all  a12 = Ai + ail  (Ai + ai2) = ail  ai2
V i = 2,..., n.
77
Hence, STEP 1 of Algorithm 3 results in a matrix M E R"' of the following form:
1
1
0
0
0
0
0
0
0
0
Since the differences azj  az , where z ranges from 1 to m, are all equal for arbitrary, but fixed j E {1,. ..,n}, all entries of M are set to 1 in STEP 2. STEP 3 yields a row vector v1 = a(1). In STEP 4, we compute a column vector u' as follows:
n
uhA (a23  J)
3=1
n
= A (a1,  aij)
j=1
= ail  all
Hence, u' = a[1]  all. STEP 5 equates v' with v1, and the algorithm finishes yielding A = u' Mv'.
Example. The following example shows that although Algorithm 3 always determines a weak decomposition of a given matrix A into outer products of column vectors and row vectors, the number of these outer products is not always minimal. In other words, there are instances when the parameter r determined by Algorithm 3 exceeds the rank of the matrix A. Let A E R4X5 be the following matrix.
1
78
2 0 0 1 4 A 5 3 1 0 8
3 5 4 2 7 3 1 0 2 6
0
= M (5 0 0 1 4)
5
4
V 0 M (5 3 1 0 8)
2
2 0
1
V 0 (3 1 0 2 0).
0
Since A is weakly decomposable into three vector pairs, the rank of rank(A) 5 3. Using Algorithm 2, we can show that rank(A) = 3. However, Algorithm 3 yields k = 4 and computes the following column vectors Ug,..., u4 and row vectors VI,..., V4 such that A = ul M V.
0 7 0 0
3 2 0 3 1 4 1
U= 5 'U2 2 ,u = 4= 3 '
12) 0) (1
V'=(2 0 2 3 2)
V2= (5 3 1 0 8)
V3 =(2 0 0 2 3),
V4 =(2 0 0 1 4) .
Computational Results. We have coded Algorithm 3 in MATLAB. We randomly generated 100 matrices of size 16 x 16 having integer values ranging from 0 to 255 and rank
79
r for r = 3,4,5. We succeeded in finding a weak decomposition of the given matrices into r vector pairs for the following number of matrices:
1. rank 3: 91%.
2. rank 4: 88%.
3. rank 5: 74%.
The average CPU time and the average amout of flops required for the rank decomposition depend on the rank of the given matrix A:
1. rank 3: 13.45 secs, 75885 flops.
2. rank 4: 18.02 secs, 110950 flops. 3. rank 5: 24.03 secs, 146710 flops. Note that the intended use of Algorithm 3 is the problem of weakly decomposing a morphological m x n template into r outer products of column templates and row templates. Savings in computation time can only be achieved if r < mn.
CHAPTER 5
THE SVD AND THE MED
5.1 Introduction
In the previous section, we discussed rankbased decompositions of matrices in the minimax algebra domain, intended for the use of morphological template decomposition. Potentially, rankbased matrix decompositions both in the linear and in the nonlinear domain could also lead to image compression techniques provided that the rank of the matrix corresponding to the image is relatively small compared to the image size. Reallife images have prohibitively large rank in linear algebra as well as in minimax algebra. However, image compression using outer product expansions does not require an exact decomposition of the corresponding matrix but merely an approximation of this matrix by a matrix of relatively small rank in terms of a suitable metric. In the linear domain, the singular value decomposition (SVD) provides the best approximation of a matrix A by a rank k matrix in the least squares sense for each k ; rank(A) [28]. Barring illconditioning of A, excellent approximations of A can be achieved even if k is significantly smaller than rank(A), making the SVD an ideal candidate for image compression purposes. Furthermore, singular value decomposition methods have proven useful in image analysis, image representation, image storage, and image restoration [1, 2].
The definition of the rank of matrices in minimax algebra in terms of outer product expansions enables us to determine a minimax algebra analogue to the SVD, which we call minimax eigenvector decomposition, or briefly MED [51]. After presenting this new
80
81
image transform, we conclude with a brief comparison of properties of the MED and the SVD.
5.2 The Singular Value Decomposition The singular value decomposition (SVD) of a realvalued m x n (m > n) matrix A allows for the representation of A as a matrix product A = UrVT, (51)
where U and V are orthogonal matrices of dimension m x m and n x n, respectively, VT denotes the transpose of V, and oi 0 0 ... 0 0 0 o2 0 ... 0 0
E= 0 0 0  0 o" (52)
0 0 0 .. 0 0
0 0 0 .  0 0
is a diagonal matrix whose diagonal entries are ordered in the form
O_ .. >O'n.
Specifically, the ai's are the square roots of the eigenvalues of the matrix ATA and are called the singular values of A. The columns of U = (u1, ... , u') are the eigenvectors of AAT and the columns of V = (vi, ... , v") are the eigenvectors of ATA. The vectors ul, ..., m are known as the left singular vectors and the vectors vi, V" are known as the right singular vectors.
A well known application of the SVD in image processing is in image compression. Observe that Eq. (51) can be rewritten as
r
A= oiui (v )T, (53)
i=1
82
where r denotes the rank of E and (vi)T denotes the transpose of the vector vi. Image compression can be achieved due to the fact that
k
A ~ Zoiui(vi)T, (54)
i=1
where k < r. Thus in order to store the whole m x n image, all we need to store is the set {0, u ,'v} .
Suppose that I IX12 symbolizes the familiar 2norm of an arbitrary realvalued matrix X. It is wellknown that lIX112 equals the largest singular value of the matrix X. The following theorem shows that a SVD of a matrix A provides the best approximation of A with respect to the metric induced by this norm.
Theorem 5.2.1. Let A = UEVt be a SVD of a matrix A e Rmxn. Setting
k
Ak = iuivi) (55)
the 2norm between A and Ak is A  AkI2 = 0k+1. Furthermore 0k+1 minimizes A B2, where B E {X E R"X : rank(X) k}.
There are many other applications of the SVD transform in image processing, pattern recognition, and linear algebra. For example, if A is of rank r < n, then the SVD provides the key for solving the least square problem of finding the vector xo that minimizes lAx  b 11. In addition, the SVD can also be used to construct the pseudoinverse of A. Employing this property of the SVD, a degraded image can be restored by multiplying the matrix modelling the degradation by its pseudo inverse [1, 2].
In spite of this great functionality, the SVD transform has not found the widespread use in computer vision as might be expected. One major reason for this is due to the complexity of the various methods for finding the eigenvalues of a matrix. The standard
83
methods for finding eigenvalues are all iterative. The amount of work required in each iteration is often prohibitively high unless initially A (or AAt) is of special form.
One of the best methods available for computing the singular values of a matrix is as follows. Observe first that if A has singular value decomposition A = UZVT and B is a matrix of form B = HAPT, where H is an orthogonal m x m matrix and P is an orthogonal n x n matrix, then B has singular value decomposition (HU)Z(PV)T. Thus the problem of finding the singular values of A can be simplified by applying orthogonal transformations to A in order to obtain a simpler matrix B with the same singular values. Golub and Kahan have shown that A can be reduced to upper bidiagonal form using Householder transformations and an excellent exposition of this method is given in [28].
This brief sketch of the methodologies for finding singular values should provide an appreciation of the complexity and intricacies of the eigenvalue problem. An additional problem concerns the accuracy of the eigenvalues one obtains by using any of the current methods. Since these iterative matrix solvers require many multiplications and divisions for large matrices, the error jui  'l between the computed and actual singular values may be quite large. In contrast, the computation of eigenvectors for minimax eigenvector decomposition does not involve these types of complexities and errors.
5.3 The Minimax Eigenvector Decomposition (MED)
An element A E R ,, is a minimax eigenvalue of A with corresponding eigenvector x E (R .)", if
A x = A +x.
Here A +x = (A + xl, ... ,A + x,), where x = (xi, ... x,).
84
Theorem 5.3.1. [14] If A E R"", then there exists a matrix U E R"' such that each column of A is an eigenvector of U corresponding to eigenvalue 0, namely
U = AN1A*.
In other words UEI A = A.
Proof. It suffices to show each of the inequalities U l A > A and U I A < A. An bound for (U M A),, can be obtained as follows:
n
(U M A), =V (uih + ahj) h=1
=y (aik +' akh) + ah
h=1 .k=1. = A (aik +' (ahk)*) + ahj h=1 .k=1. >A(aik +' (aik)*) +aij k=1
> ai
Vi=1,...,m;Vj=1,...,n. The last inequality holds because for any a E R ,, the sum a +' (a)* equals 0 for finite ro V a ER
a and oc else. Similarly, a + a* = a* + a = oo else . This fact and the first associative inequality in Theorem 3.2.2 lead to a proof of U M A < A similar to the
85
above proof of the inequality U M A > A.
(U M A)j =((A 9 A*) M A)i,
(A (A* M A)),J n n I
= aih (ahk + akj)
h=1 k=1
n n
= aih +1 V(akh)* + akj0l
h=1 .k=1
n
<5 aij ' ((akj)* + ak) k=1
aij
V I = 1, ..., m; V J = 1,...,n.
Therefore, U M A = A and we have certain parallels to the SVD of linear algebra. Specifically, the MED of A has form
A = U M f M Vi,
(56)
where V = At,
0*1
00
 00
00
00 and a = 02 = =on =0.
00
02
00
00
00
00
00
00
00
00
Equivalently,
A = (A Z A*) M A.
It follows from Eq. (56) that there are hardly any computations involved in finding the MED of a matrix. The eigenvalues are already known to be zero. One only needs
00
00
00
00
00
00
00
On
00
00
86
to form the conjugate matrix A* and compute the product U = A V9 A*. In analogy to Eq. (53), it follows from Eq. (56) that
A = 10i + (Ui M (vi),) . (57)
i=1
Since (v')t denotes the ith row vector a' of A and oi = 0 for i = 1,2, ... ,n, we actually have the simpler formulation
n
A = (Ui M ai) . (58)
1=1
Since the ui's are the columns of the matrix U = A E A*, the determination of the components of the ui's requires only additions and the logic operation A of minimization. The components of the row vectors a of A are given and the computation of the eigenimage u' M a' requires only additions. Thus the determination of eigenvectors and eigenimages required in the minimax decomposition of a matrix is for all practical purposes trivial. The following example serves to demonstrate this observation.
Suppose
= 3 7 8 3 9 5)
A=(9 31 he * 78 ( 4 I)
5 1 0 8 7 0
and, therefore,
15 8)
U = AIAJA* =(60 1 .8
(2 7 0
5.4 Comparison of the MED and the SVD
In comparison to the SVD, minimax decomposition of a matrix is trivial. However, what is often desired is not the complete set of outer product expansions but only a small
87
subset of these with the property that this subset provides for the best approximation Ak of the original matrix A. In case of the SVD, once the singular values and singular vectors for the series expansion of Eq. (53) have been determined (the difficult part), it is trivial to obtain the best approximating subset of k outer product expansions for a given k < n with respect to the 2norm since by Theorem 5 k k
A  ori iu (Vi)T < A  Ci cii(yi)T
2 i=1 2
Vci E R, xi E R, Vyi E R (i=1,...,k). In case of the MED, the determination of the eigenvectors for the maxseries expansion of Eqs. (57) or (58) is trivial, and the eigenvalues oi are known to be 0. However, even the selection of the best approximating subset of k eigenvector pairs with respect to a given metric d is a difficult problem. This problem consists of finding a set of indices I {1,..., n} with cardinality k such that d (A,VU.M (vi) d (A, V uM(v') T)
\iEI / \iEP'
V' P {C ..n : I' = k . Clearly, the straightforward approach requires ( comparisons which are prohibitively many. Thus, we decided to use the following iterative algorithm to select a set of k eigenvector pairs. We chose d, as a metric in order to avoid multiplications. Other convenient choices of the distance d would be d2 and the metric induced by the Frobenius norm.
1. Find two vectors uj and aj such that
n
d, (A, u'. M a) A di (A, u' M a') i=1
Rename this pair as ul and al; i.e. permute 1 and j
88
2. Suppose the first j (J > 1) eigenvector pairs have been found. Find two vectors
ua and a' such that
d, (Ua M a"r) V V (u' M a') , (u' a')
n
di( (u M a') V \ (u M a') ,V (u M al)
i=+l==
Rename this pair as uj+1 and aj+1.
3. Set j = j + 1 and continue with step (2.) and step (3.) until j = k.
k
An approximation of A with respect to the norm di is given by Ak = V (u' 1=1
M a'). Note
that if A E Zmxn, then all necessary computations can be performed within Z. Hence, this algorithm guarantees numerical stability for the case A E Z"'.
Example. If we apply Algorithm 4 to the 5 x 5 matrix /39 7 0 94 84\ 62 0 58 11 3 A= 73 19 30 2 01,
9 75 8 20 29 \58 27 18 82 4/ then using the first three eigenvector pairs we obtain the reconstructed matrix
3
A3= V (uM a')=
i=1
39 7 62 0 34 0 9 0 \58 27
0 58 30
5
18
94 11
2
20 82
84
3
0 10
4 /
Using the complete set of five eigenvector pairs we obtain, of course, A5 = A.
3T
To contrast this with the SVD, the series expansion A3 = OI aui(vi) yields i=1
/36.813789
69.110519 A3 = 63.394817
10.670452
\60.955730
8.839371 1.802377
19.603844 75.780029 23.081055
2.534787
43.164482 41.031982 3.457618 25.672878
103.915436 15.181345 10.571262 22.650938 61.892658
70.123454
7.691323
7.335522 23.932810 33.012875 /
89
while
/39.000019 6.999990 0.000003 94.000031 84.000000
62.000015 0.000008 58.000008 11.000007 2.999976
A5 = 73.000015 18.999996 30.000013 1.999997 0.000021
9.000000 75.000031 8.000004 20.000019 29.000017
\58.000008 26.999996 18.000004 82.000008 4.000021 /
Note that although in theory A5 = A for the SVD as well, in actual computation A5 is only approximately equal to A due to numerical errors introduced by finite digital computers when using multiplications and divisions.
Example. As an additional example, we compare the MED and SVD transform on an actual image. Figure 9 shows the 12 x 12 input image, while Figures 10 and 11 show the eigenvector expansions for three and seven eigenvector pairs of the SVD and the MED, respectively.
Figure 9 The 12 x 12 source image A.
Figure 10 Image reconstruction using the SVD transform; the image on the left corresponds to A3 and the image on the right to A7.
90
Figure 11. Image reconstruction using the MED transform; the image
on the left corresponds to A3 and the image on the right to A7.
Examples like these indicate that image reconstruction using the MED works well for "boxy" images of relatively small size. However, if A E Rmx" corresponds to an m x n real life image, the distances d, (A, Ak) are prohibitively large for all k < n. Moreover, Corollary 4.6.2 implies that the matrix Ak does not necessarily represent the best approximation with respect to the metric d, of A among all matrices of rank k.
We dealt with these problems concerning the image reconstruction of reallife images using the MED by employing the technique of divide and conquer. Let A E Rm'xn. If p divides m and q divides n, then A has the following block structure, where A' E Ri x : A1 A2 ... A
Aq+1 Aq+2 2q
A= . . A)
A ipq+1 . ..Apq
For given k < n and for each i = 1,. . . , pq, an application of Algorithm 4 to A' E R P X yields a matrix (A'), E R x i, which we conveniently denote by A'. Let A' A2 ... Aq
A(pq) _ A+ A +2 ... A2q
AI . .+1 . A.1
(A(1 q~i ... ... A /2
pq
Note that d, (A, Ak) = di (A, Ak).
i=1
91
Example. Let us apply this technique to the following realvalued image of size 32 x 32, which corresponds to a matrix A E R32X32. The matrix A' can be stored on a digital computer using only a quarter of the entries of the matrix A E R32x32, since the representation of the matrix A' = uM[O (vi)T E R8x8 requires only 2  8 =16 entries (4,4)
compared to the 64 entries of an 8 x 8 matrix. Similarly, forming A,' yields an image compression rate of 50 %, since the matrices A' E RS 8 have representation in terms of
2 outer products of vector pairs of length 8.
Figure 12. The 128 x 128 source image A.
N
Figure 13. Image reconstruction using the SVD transform; the image on the left corresponds to A16 and the image on the right to A32.
92
Figure 14. Image reconstruction using the MED transform; the image on the left
corresponds to about one fourth and the image on the right to half of the total information.
CHAPTER 6
MORPHOLOGICAL DECOMPOSITION AND APPROXIMATION
6.1 Introduction
The previous chapter describes the MED transform and it's use for determining the best approximation of a matrix with respect to the metric di. By Corollary 4.6.2, the outer product expansions given by the MED transform do not necessarily yield the best approximations of the original matrix.
In this chapter, we are presenting an integer programming approach to matrix approximation problems using d, as a distance. By the definition of a metric, matrix approximation in terms of outer product expansion results in an exact representation of the original matrix, if the distance between the norm and the approximating outer product expansion is 0. Due to the natural correspondences between matrices and images as well as between matrices and templates, solving matrix approximation problems could have dual use for image and template applications. However, since the size of images would lead to immense optimization tasks, this approach is only intended for the purpose of template decomposition and approximation. Using this approach, we are able to tackle the problem of weakly decomposing arbitrary given templates with respect to the operation M in its most general form. We provide exact definitions below.
6.2 Morphological Decomposition and Approximation Problems
The decomposition problems, we are primarily interested in, are matrix and template decompositions with respect to the operation M . Furthermore, we assume that the
93
94
templates which are to be factorized are elements of (Rx,)X where X is a finite subset of Z2 (without loss of generality X contains z := (0, 0). Linear mixed integer programming is suited to solve a very general form of these problems. The General Form of the Decomposition Problem. Let A E Rmx". Assume that A is decomposable as follows for some natural number k: A = B M C,
B e Rmxk
C kxn.
C E Rix
Suppose a set Y is defined as follows: Y = {Yi, ... ,Ym*k+k*n} = {b1,. .., blk, ..., bmi,... bmk, C11, ..., CI,.. CkI, ..., Ckn}
Then there are linear functions flj (i = 1, .. ., m; j = 1,... ,n; 1 = 1,... ,k) such that for each value of i and j the following inequalities hold: fij(yi, ...,Ymk+kn) < aij f 2;j(Y1, ... , ymk+kn) < aij
(ij)
f(yi,. . ., Ymk+kn) < aiMoreover, for each i, j there exists an index l" E {1,..., k} such that: fs(yi, . . ., Ymk+kn) aij
Obviously the linear functions f! are simply derived from the additive maximum B M C:
f: I(yI, . . ., Ymk+kn) = bil + cj.

Full Text 
PAGE 1
MATRIX DECOMPOSITION IN MINIMAX ALGEBRA AND APPLICATIONS IN IMAGE PROCESSING By PETER SUSSNER A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1996 UNIVERSITY OF FLORIDA LIBRARIES
PAGE 2
ACKNOWLEDGMENTS First and foremost, I would like to thank Dr. Gerhard X. Ritter for providing me with excellent opportunities to perform interdisciplinary research on his contract and for guiding my learning process as a mathematician. I am grateful to Dr. Panos M. Pardalos for giving new direction and inspiration to my research. I thank all my committee members, in particular Dr. David C. Wilson, for many stimulating and helpful conversations. I acknowledge Dr. Hongchi Shi, Dr. Joseph Wilson, Mrs. Monica Sweat, and all the faculty and graduate students at the Departemt of CISE who assisted me in many technical and computer science related questions. My father, Walter Sussner, deserves my deepest gratitude for his neverending support and encouragement. I am also indebted to Dr. Patrick Coffield and Dr. Sam Lambert from the Air Force Armament Laboratory for funding this research under US Air Force Contract F0863589C0134. Finally, I acknowledge the American and the German taxpayer for their investment in my education. I hope that I will be able to justify their financial support through my future contributions to society. ii
PAGE 3
TABLE OF CONTENTS ACKNOWLEDGMENTS ii ABSTRACT v CHAPTERS 1 INTRODUCTION 1 LI Motivation 1 1 .2 Introduction to Template Decomposition 2 1.3 Summary of Results and Organization 6 1.4 Notations 7 2 IMAGE ALGEBRA BACKGROUND 9 2.1 Introduction 9 2.2 Operands 10 2.3 Operations on Images 12 2.4 Operations between Images and Templates 13 2.5 Operations between Templates 15 2.6 Metrics on Matrices and Invariant Templates with Finite Support 17 2.7 Decompositions of Templates 18 2.8 Relevance of Template Decompositions 18 3 MINIMAX ALGEBRA CONCEPTS AND EXTENSIONS 20 3.1 Introduction 20 3.2 Belts and Blogs 21 3.3 Spaces of Vectors and Matrices 26 3.4 Conjugacy 30 3.5 The Principles of Opening and Closing 32 3.6 Linear Independence and Rank 35 4 RANKBASED DECOMPOSITIONS OF MORPHOLOGICAL TEMPLATES . 45 4.1 Introduction 45 4.2 Correspondence between Matrices and Rectangular Templates 46 4.3 Invariance Properties of the Rank of Realvalued Matrices 49 4.4 Separable Matrices and Templates 51 4.5 Matrices and Templates having Rank 2 56 4.6 Matrices and Templates of Arbitrary Rank 64 iii
PAGE 4
5 THE SVD AND THE MED 80 5.1 Introduction 80 5.2 The Singular Value Decomposition 81 5.3 The Minimax Eigenvector Decomposition (MED) 83 5.4 Comparison of the MED and the SVD 86 6 MORPHOLOGICAL DECOMPOSITION AND APPROXIMATION 93 6.1 Introduction 93 6.2 Morphological Decomposition and Approximation Problems 93 6.3 Integer Programming Formulations 97 7 COMPLEXITY OF MORPHOLOGICAL DECOMPOSITION PROBLEMS . 105 7.1 Introduction 105 7.2 An NPComplete Problem 105 7.3 Extensions 110 8 CONCLUSION AND SUGGESTIONS FOR FURTHER RESEARCH .... 112 REFERENCES 115 BIOGRAPHICAL SKETCH 120 iv
PAGE 5
Abstract of Dissertation Presented to the Graduate School of the University of Rorida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy MATRD( DECOMPOSITION IN MINIMAX ALGEBRA AND APPLICATIONS IN IMAGE PROCESSING By Peter Sussner August 1996 Chairman: Dr. Gerhard X. Ritter Major Department: Mathematics Minimax algebra is a mathematical theory which originated from problems in operations research and machine scheduling. Surprising parallelisms exist between linear algebra and minimax algebra. Both linear algebra and minimax algebra can be embedded into the mathematical theory of image algebra which forms a unified mathematical environment for all image processing and computer vision applications. Methods for matrix decomposition in linear algebra have found numerous applications in image processing, in particular for the problem of template decomposition. Therefore, it seems reasonable to investigate matrix decomposition techniques in minimax algebra with applications in image processing. We provide a mathematical basis for these investigations by establishing a new theory of rank within minimax algebra. Thus far, V
PAGE 6
only minimax decompositions of rank 1 matrices into outer product expansions are known to the image processing community. This dissertation derives various new methods for the decomposition of matrices of rank larger than 1 within minimax algebra. Moreover, this thesis proves the NPcompleteness of the problem of decomposing arbitrary matrices in minimax algebra. Consequently, we obtain a fundamental result for the area of morphological image processing: the NPcompleteness of the general morphological template decomposition problem. vi
PAGE 7
CHAPTER 1 INTRODUCTION Â» . . . I'1.1 Motivation A number of researchers have discovered that a variety of problems in operations research and mathematical economy, e.g. some problems in optimization on graphs and networks, in machinescheduling, and approximation theory, can be formulated using a certain algebraic lattice structure [3, 5, 8, 13, 27]. We can think of this algebraic structure as the real numbers (perhaps extended by the symbols +oo and/or Â— oo) together with operations maximum and/or minimum. In the problem formulations, which we referred to above, similarities to concepts from linear algebra become apparent. In analogy to linear algebra, CuninghameGreen presented the minimax algebra, a unified account of the algebra of spaces of vectors over a lattice [14]. Jennifer Davidson showed how to embed the minimax algebra into the image algebra [15, 16]. The theory of image algebra incorporates all known image processing techniques in a translucent manner [19, 21, 49, 47, 52]. One powerful implication of this is that all the tools of minimax algebra are directly applicable to solving problems in image processing whenever the minimax product of matrices, denoted by El , is involved. Jennifer Davidson's embedding of the minimax algebra into the image algebra maps vectors to images and it maps matrices to templates. This dissertation makes use of two different mappings. One of these mappings relates matrices to images, and the other one relates matrices to rectangular templates. One particular consequence of this relationship between matrices and images is that the Singular Value Decomposition known from 1
PAGE 8
numerical linear algebra can be used for the purpose of image restoration and image compression [1,2, 35]. On the other hand, matrix decomposition methods lead directly to methods for template decomposition. The problem of template decomposition Â— both in the linear and in the nonlinear domain Â— has been extensively investigated by a host of researchers. We present a brief review of different approaches to this problem below. A multitude of researchers have applied matrix decomposition techniques stemming from linear algebra to image processing. However, the minimax algebra domain has been widely neglected in this respect. Since operations derived from minimax algebra also play an important role in image processing, we consider it timely to present an exposition on the theory of matrix decomposition in minimax algebra with applications in image processing. 1 .2 Introduction to Template Decomposition The theory of mathematical morphology grew out of the belief that many image processing operations can be expressed in terms of a few elementary operations [32, 33, 41, 54, 55]. The origins of mathematical morphology lie in the work done by H. Minkowski and H. Hadwiger on geometric measure theory and integral geometry [31, 42, 43]. Both linear convolution and morphological methods are widely used in image processing. One of the common characteristics among them is that they both require applying a template to a given image, pixel by pixel, to yield a new image. In the case of convolution, the template is usually called convolution window or mask; while in mathematical morphology, it is referred to as structuring element. Templates used in realizing linear convolutions are often referred to as linear templates. Templates can vary greatly in their weights, sizes and shapes, depending on the specific applications.
PAGE 9
3 Intuitively, the problem of template decomposition is that given a template t, find a sequence of smaller templates t^, . . . ,tÂ„ such that applying t to an image is equivalent to applying tj, . . . , tÂ„ sequentially to the image. In other words, t can be algebraically expressed in terms of ti, . . . ,tn. One purpose of template decomposition is to fit the support of the template (i.e. the convolution kernel) optimally into an existing machine constrained by its hardware configuration. For example, ERIM's CytoComputer [56] cannot deal with templates of size larger than 3 x 3 on each pipeline stage. Thus, a large template, intended for image processing on a CytoComputer, has to be decomposed into a sequence of 3 x 3 or smaller templates. A more important motivation for template decomposition is to speed up template operations. For large convolution masks, the computation cost resulting from implementation can be prohibitive. However, in many instances, this cost can be significantly reduced by decomposing the masks or templates into a sequence of smaller templates. For instance, the linear convolution of an image with a grayvalued n x n template requires n? multiplications and 1 additions to compute a new image pixel value; while the same convolution computed with an 1 x n row template followed by an n x 1 column template takes only 2n multiplications and 2(n 1) additions for each new image pixel value. This cost saving may still hold for parallel architectures such as mesh connected array processors [36], where the cost is proportional to the size of the template. The problem of decomposing morphological templates has been investigated by a host of researchers. Zhuang and Haralick [62] gave a heuristic algorithm based on tree search that can find an optimal twopoint decomposition of a morphological template if such a decomposition exits. A twopoint decomposition consists of a sequence of templates each consisting of at most two points. A twopoint decomposition may be best suited
PAGE 10
4 for parallel architectures with a limited number of local connections since each twopoint template can be applied to an entire image in a multiplyshiftaccumulate cycle [36]. Xu [60] has developed an algorithm, using chain code information, for the decomposition of convex morphological templates for twopoint system configurations. Again using chaincode information, Park and Chin [46] provide an optimal decomposition of convex morphological templates for 4connected meshes. However, all the above decomposition methods work only on binary morphological templates and do not extend to grayscale morphological templates. A very successful general theory for the decomposition of templates, in both the linear and morphological domain, evolved from the theory of image algebra [19, 21, 47, 52, 48] which provides an algebraic foundation for image processing and computer vision tasks [49]. In this setting, Ritter and Gader [21, 50] presented efficient methods for decomposing discrete Fourier transform templates. Zhu and Ritter [61] employ the general matrix product to provide novel computational methods for computing the fast Fourier transform, the fast Walsh transform, the generalized fast Walsh transform, as well as a fast wavelet transform. In image algebra, template decomposition problems, for both linear and morphological template operations, can be reformulated in terms of corresponding matrix or polynomial factorization. Manseur and Wilson [39] used matrix as well as polynomial factorization techniques to decompose twodimensional linear templates of size m x n into sums and products of 3 x 3 templates. Li [37] was the first to investigate polynomial factorization methods for morphological templates. He provides a uniform representation of morphological templates in terms of polynomials, thus reducing the problem of decomposing a morphological template to the problem of factoring the corresponding polynomials. His approach provides for the decomposition of onedimensional morpho
PAGE 11
5 logical templates into factors of twopoint templates. Crosby [11] extends Li's method to twodimensional morphological templates. Davidson [18] proved that any morphological template has a weak local decomposition for meshconnected array processors. Davidson's existence theorem provides a theoretical foundation for morphological template decomposition, yet the algorithm conceived in its constructive proof is not very efficient. Takriti and Gader formulate the general problem of morphological template decomposition as optimization problems [23, 58]. Sussner, Pardalos and Ritter [57] use a similar approach to solve the even more general problem of morphological template approximation. However, since these problems are inherently NPcomplete, researchers try to exploit the special structure of certain morphological templates in order to find decomposition algorithms. For example, Li and Ritter [38] provide very simple matrix techniques for decomposing binary as well as grayscale linear and morphological convex templates. A separable template is a template that can be expressed in terms of two onedimensional templates consisting of a row and a column template. Gader [22] uses matrix methods for decomposing any grayscale morphological template into a sum of a separable template and a totally nonseparable template. K the original template is separable, then Gader' s decomposition yields a separable decomposition. If the original template is not separable, then his method yields the closest separable template to the original in the mean square sense. Separable templates are particularly easy to decompose and the decomposition of separable templates into a product of vertical and horizontal strip templates can be used as a first step for the decomposition into a form which matches the neighborhood configuration of a particular parallel architecture. In the linear case, separable templates are also called rank 1 templates since their corresponding matrices are rank 1 matrices. O'Leary [44] showed that any linear template of rank 1 can be factored exactly into a
PAGE 12
product of 3 X 3 linear templates. Templates of higher rank are usually not as efficiently decomposable. However, the rank of a template determines upper bounds of worstcase scenarios. For example, a linear template of rank 2 always decomposes into a sum of two separable templates. 1.3 Summary of Results and Organization In nonlinear or morphological image processing, convolutions of an image with a template consist of combinations of maximum, minimum, and addition operations. The theory of minimax algebra provides the mathematical background for the algebraic structures arising from these operations. This dissertation presents a new concept of matrix rank within minimax algebra which is similar to the one in linear algebra. Due to the correspondence of matrices and rectangular templates, we obtain a concept of morphological template rank. In this context, separable morphological templates represent exactly the templates of rank 1. A separable template is easy to decompose into an outer product of a column and a row template. The cost savings in computational complexity achieved by a decomposition of this form still remain, if templates of arbitrary, but fixed rank k are decomposed into outer products of column and row templates. First, we present an algorithm which, for any rectangular morphological template of rank 2, constructs a weak decomposition of this template into two column templates and two row templates in polynomial time. Afterwards, we present a heuristic algorithm for the decomposition of arbitrary templates of relatively small rank k. In linear algebra, an exact decomposition of a rank k matrix into k matrices of rank 1 can be determined by using the SVD. This thesis provides a matrix decomposition in minimax algebra which is mathematically very similar to the SVD in linear algebra. However, this decomposition which we call MED, generally does not provide the best
PAGE 13
7 approximation of a matrix by A; outer products of column and row templates in terms of an appropriate matrix norm. Instead, we solve this problem by formulating it as an integer programming problem. Furthermore, the general morphological template approximation problem can be formulated as an integer programming problem. We provide preliminary computational results. Finally, we justify these integer programming approaches by showing that the general morphological template decomposition problem is NPcomplete. Indeed, the problem of weakly decomposing arbitrary rectangular templates into column and row templates is already NPcomplete. This dissertation is organized as follows: First, we introduce the reader to the language of image algebra. This mathematical theory is suited to describe all image processing operands and operations in a translucent manner. Since we focus on morphological image processing, we proceed with a brief description of the algebraic structures defined in minimax algebra. We relate matrices in minimax algebra to rectangular morphological templates. In chapter 3, we establish a rank method for the morphological decomposition of matrices and rectangular templates. We briefly review the Singular Value Decomposition in chapter 4. We derive a similar decomposition in minimax algebra, called Minimax Eigenvector Decomposition. In chapter 5, new integer programing approaches for solving arbitrary instances of the general morphological template decomposition problem are presented. Finally, we justify this strategy by showing that the problem of decomposing arbitrary morphological templates is NPcomplete. 1.4 Notations Before we proceed, let us discuss some frequently used notations. As usual, the symbol N stands for the set of natural numbers, Z for the set of integers, and R stands
PAGE 14
8 for the set of real numbers. The symbol R"*" denotes the positive real numbers. The real numbers extended by the symbol Â— oo are denoted by R_oo and the real numbers extended by the symbol +00 are denoted by R+ooSimilarly, Rl^^ stands for R"*" U {00}, and R'\.oo stands for R"*" U {+00}. When adjoining both +00 and Â—00 to R, one obtains the set RÂ±ooWe use the symbols V and A to denote the binary operations of maximum and minimum. The positive part of a real number x, denoted by is defined as 0 V x. Similarly the negative part of a: is = 0 V {Â—x). Matrix inequalities are always defined entrywise. For example, given matrices A and B, the inequaHty A < B means aij < bij If 7 is an associative and commutative operation on a finite set X = {xi, . . . ,xÂ„}, then we abbreviate the expression 0:17 . . . 'yxn by r"_j(a;,) or by.F^ ^a:,. The symbol Tx is yet another way to denote 0:17 . ..^Xn provided that x denotes the transpose of (xi,...,xÂ„), or briefly (xi, . . . ,a;Â„) . If / is a function on X, then f{X) denotes {/(xi),...,/(xÂ„)}. A sequence A;2, Â• Â• Â• , A;Â„) is called strictly increasing if ki < k2 < ... < K.
PAGE 15
CHAPTER 2 IMAGE ALGEBRA BACKGROUND , 2.1 Introduction Images not only provide important sources of information in every day life, but also in various sciences and in the military. However, often large parts of the information contained in an image remains hidden from the eye of the human observer. Image processing techniques, when implemented on a computer, can help to filter out these parts of the information. For instance, image processing methods can be applied in order to identify and locate objects in an image. The objects the user is looking for can be as diverse as cancer cells in medical applications and tanks in military applications. In cardiology, for example, physicians need to find the location of the anterior and posterior wall of the heart. A multitude of image processing methods for image smoothing, image enhancement, edge detection, thinning, skeletonizing etc. exist [4, 30, 29]. When faced with a specific problem, like finding certain objects in an image, the developer of an image processing algorithm has to make the difficult decision which techniques to use. In this case, a certain amount of a priori knowledge, e.g. on the shape and size of the objects, will turn out to be very useful. For a long time developers of image processing software have documented their algorithms by using lengthy word descriptions and analogies that are extremely cumbersome and often ambiguous. These ad hoc approaches led to a proliferation of nonstandard notation and increased research and development cost. In response to this chaotic situation, 9
PAGE 16
10 G. X. Ritter et al. have developed the (AFATL) image algebra since the mid1 980' s. This algebraic structure is specifically designed for image manipulation. The image algebra constitutes a heterogeneous or manyvalued algebra in the sense of Birkhoff and Lipson [7, 48]. Intuitively, an algebra is simply a collection of nonempty sets together with finite number of mappings (operations) from one of the sets into another. In image algebra, six basic types of operands exist, namely value sets, coordinate sets, the elements of each of these sets, images and templates. In this dissertation, we restrict our attention to singlevalued images. We refer the reader to other publications for a discussion of multivalued images [48, 52]. 2.2 Operands A point set is a subset of a topological space. In digital image processing, point sets are often assumed to be rectangular discrete arrays. A value set is a homogeneous algebra such as Z, R, Roo(= R U {oo}), R+oo(= R U {oo}) or C. An Fvalued image a is the graph of a function a : X Â— > F where X denotes some point set and F denotes some value set. Obviously, a can be represented as a matrix if X is a rectangular discrete array. The set of all images a : X ^ F is denoted by F^. The notation R^oo is frequently used instead of (R_oo)^. Let X, Y be point sets and F a value set with neutral element 0. A generalized Fvalued template t from Y to X is a function t : Y Â— > F'^. The set of all templates t : Y Â— > F'^ is denoted by (F'^)^. So for each y G Y the value t(y) is an Fvalued image on X. We define ty = t(y) and call the point y the target point of t. The function values ty(x), where x G X, are called weights of the template t at location y.
PAGE 17
11 The support of ty is defined as follows: 5(ty) = {xeX : ty(x)^0} Note that the neutral element 0 depends on the algebraic structure of F. For example, the element 0 represents the neutral element for R together with +, whereas the element oo represents the neutral element for R_oo together with V. A template t 6 (p^) is called translation invariant if and only if for each triple X, y, z E X with x + z G X and y + z G X the following equation holds: ty(x) = ty+z(x + z) . In this case, it suffices to know ty^ for an arbitrary yo Â€ X. If X is discrete and 5(tyo ) is finite then t can be represented pictorially. For example, let X = Z'. Moreover, let r 6 (R^oo) ^ the translation invariant template which is determined at each point y = (a;,?/) G X by the following function values of x e X: i/ x = y if x = y(0, 1) z/ x = y + (0, 1) ifx\Â—x Â— \ and y Â— l
PAGE 18
12 y\ Â— 1 Â— y+1 Â— 1 Â— Â— 1 Â— JC1 2 5 10 X 4 / \ nV 12 JC+1 3 0 5 Figure 1. The support of the template r at point y. The hashed cell indicates the location of the target point y = (a;,t/). A translation invariant template r is called rectangular, if 5(ry) forms a rectangular discrete array. We speak in particular of an m x n template, if 5'(ry) has size m x n. An m X 1 template is also called a column template of length m, and a 1 x n template is called a row template of length n. For the remainder of this dissertation, point sets can always be identified by the symbols X, Y, Z and value sets by the symbols F,G, H. The point set X is always assumed to be finite. 2.3 Operations on Images The basic operations on Fva!ued images stem from the naturally induced operations on the value set F. In particular, addition, multiplication, and maximum are defined as follows: a + b = {(x, c(x)) : c(x) = a(x) + b(x), Vx G X} V a, b G R^^o a Â• b = {(x, c(x)) : c(x) = a(x) Â• b(x), Vx Â€ X} V a, b G a V b = {(x, c(x)) : c(x) = a(x)Vb(x), Vx Â€ X} V a, b Â€ R^oo Â• We adopt the usual convention r + (Â— oo) = (Â— oo) + r = Â— oo for all r e Roo
PAGE 19
13 2.4 Operations between Images and Templates The more complex operations which we define in this section are based on a binary operation o : F x G ^ H and an associative and commutative operation 7 on the value set H. If a e G^ and t G (F^) , then a@t 6 Â— the forwards generalized product of the image a with the template t Â— can be formed: , aÂ®t={(y,b(y)) : b(y) = F (a(x) o My)), y Â€ y} . Hence, a forwards product of this form can effect a transformation from one image having a certain value set and a certain point to another image having a different value set and a different point set. This thesis concentrates on the backwards products of images and templates which are defined below. If a 6 F^ and t G (G^)^, then a@t e Â— the (backwards) generalized product of the image a with the template t Â— is defined as the following image: aÂ®t = {(y, b(y)) : b(y) = r^(a(x) o ty(x)), y Â€ y} . Specific examples of this operation are given below. Unless stated otherwise, we are referring to backwards imagetemplate products for the remainder of this thesis. If, for example, a 6 R^oo * ^ (f^?oo)^> then the generalized convolution ("0") and the additive maximum (" El ") are defined as follows: a0t= i(y,b(y)) : b(y) = \J a(x) + ty(x), y Â€ Y I xex a et = I (y, b(y)) : b(y) = ^ a(x) * ty(x), y Â£ Y V xex
PAGE 20
14 where \J a(x) + ty(x) = max{a(x) + ty(x), x G X} . xex The operation additive minimum (" El ") is defined in a similar fashion. Many image processing techniques such as the Rolling Ball Algorithm and algorithms for noise removal employ imagetemplate products of the form aElt or aOt. Both forms appear in the MaxMin Sharpening Algorithm, which enhances blurred images. 1 ft m Figure 2. Blurred characters and result of applying maxmin sharpening to blurred characters One of the most widely used transforms in image processing is the Fourier transform [10, 59]. Like many other image processing transformations, the Fourier transform can be expressed in terms of a generalized convolution. This means that there is a template f such that the Fourier transform of an image a is given by a 0f . The images in Figure 4 show an image of a jet and its corresponding Fourier transform image after centering. Notice that almost all of the information is contained in the center of the latter image.
PAGE 21
15 1 Figure 3. The image of a jet and its centered Fourier transform image. Typical applications of the Fourier transform are low and high pass filtering [29, 49]. The figures below demonstrate how low pass filtering can be applied to medical images. Figure 4. Noisy angiogram and filtered angiogram using a low pass filter 2.5 Operations between Templates The basic operations between generalized templates are induced by the basic operations on images. Thus, the addition, the multiplication, and the maximum of two templates are straightforward: s + t = {(y, ry) : ry = Sy + ty, y G Y} Vs,tG (r^^) s*t = {(y, ry) : ry = Sy * ty, y Â€ Y} Vs,te (r^) Y sVt = {(y,ry) : ry = Sy V ty, y 6 Y} Vs,tG (r?^)
PAGE 22
16 Since the set R_oo together with the morphological operations and V forms an algebraic structure in minimax algebra, the R_oo valued templates are also called morphological templates. Templates can not only be convolved with images, but also with other templates. Ritter [48, 52] defines the generalized product of a template s e (F^) and a template t G (G^) , also denoted by s Â®t, as the template r e (H^) by setting ry(z) = r^(sx(z) o ty(x)) Vy G Y, Vz e Z . For example, realvalued templates t and s allow us to form a generalized convolution ("@") and an additive maximum ("0"): (s et)y(z) = * ty (^)) Vy G Y, Vz Â€ Z xGX (sElt)y(z)= \/{sM + hi^)) VyeY, VzeZ XGX The following column templates r, s, t G (Roo) satisfy r = sEI t. r = y Sy = 0 , Figure 5. The template r constitutes the additive maximum of the templates s and the template t. Note that the generalized products s Â®t and t @s of an m x 1 rectangular template s and a 1 X n rectangular template t always lead to m x n rectangular templates.
PAGE 23
17 Suppose that X = and that s, t e (R^) are the following translation invariant templates: The template product r = s 0 1 is given by the template r in Figure 1. 2 tv = 7 12 \ / Figure 6. Pictorial representation of a column template s and a row template t. 2.6 Metrics on Matrices and Invariant Templates with Finite Support It is well known that a metric on a metric space M induces a metric on the space jVf""^" of all m X n matrices. Let d be the metric on R given by: d{a, b) = \a b\ Va, 6gR. Then the following definitions induce metrics on R"Â»><": n m i=\ j=l d2iA, B) = \l d{a,j, hij) . hi Similarly, if s Â€ (R''^)^ and t G (R^)^ are invariant templates with the same finite support S at some target point y E Y, then the following equations provide metrics on the set of all invariant templates with that particular finite support: ^l(S, t) = ^4ty(x), Sy(x)) xes C'2(S, t)= y d{ty{x), Sy(x)) . xeS Since the templates s and t are invariant, these definitions do not depend on the choice of the target point y Â€ Y.
PAGE 24
18 2.7 Decompositions of Templates A sequence of templates (t\ . . . , t^") in (F^) together with a strictly increasing sequence of natural numbers (fci, . . . , /;Â„) is called a weak decomposition (with respect to the operation "Â®") of a template t 6 (F^) if the template t can be represented as follows: t = (t^ Â® . . . (t^^+^ Â® Â• Â• Â®t*=^)7 Â• Â• Â• 7 (t'="^+' Â® . . . Â®t^") . In the case where n = 1, we speak of a strong decomposition of the template t. For example, a strong morphological decomposition of a template t E (R?oo) has the following form: t = 0 e 0 . . . sa 1^= . If the template t strongly decomposes into a column template and a row template t^, it is referred to as a separable template. For example, the template r E (R?oo) given in Figure 1 represents a separable template since this template decomposes into the column template s 6 (R?oo) ^"d the row template t G (R?oo) of Figure 6. 2.8 Relevance of Template Decompositions The following associative and distributive laws hold for an arbitrary image a 6 F^ and arbitrary templates t e F'^ and s e F'^: aEl (sEI t) = (aEI s) Elt , a0 (s Vt) = (aEls) V (aElt) .
PAGE 25
19 These results establish the importance of template decomposition. Realizing the convolution of an image a with a template t on a computer is an extremely complex task, since each pixel value at location y depends on all pixel values of a within the neighborhood given by ty. Decomposing the template t often helps to reduce the number of computations involved in forming an imagetemplate product by providing alternative ways to compute this convolution. For example, the application of a separable template to an image is equivalent to the sequential applications of a column template and a row template. Hence, a quadratic problem reduces to a linear problem. Computing the Fourier Transform of an image represents a linear convolution. The speedup of this computation achieved by the Fast Fourier Transform relies heavily on linear template decomposition. Template decompositions have also been successfully used for the purpose of achieving a better fit to the computer architecture. For example, if each processing element has direct neighbors in N, S, W, and E directions, then the user of image processing software benefits from template decompositions into invariant templates having supports of the following form: / \ \ /
PAGE 26
CHAPTER 3 MDSflMAX ALGEBRA CONCEPTS AND EXTENSIONS 3.1 Introduction In analogy to linear algebra, CuninghameGreen presented the minimax algebra, a unified theory of the algebra of spaces of vectors over F [14]. Minimax algebra includes familiar concepts from linear algebra such as linear dependence of vectors, matrix rank, and eigenvalues of matrices. Operations on vectors and matrices in minimax algebra resemble the operations on vectors and matrices in linear algebra. In his Ph.D. dissertation, Paul Gader embedded the linear algebra into the image algebra [20]. Likewise, J. L. Davidson provided an isomorphism of the minimax algebra into the image algebra [15, 17]. This isomorphism maps vectors to images and matrices to templates. Hence, results on matrix decomposition in minimax algebra induce results on template decomposition in image algebra. This dissertation exploits a different kind of relationship between matrices in minimax algebra and rectangular morphological templates, which is also useful for the purpose of template decomposition. A similar relationship exists between matrices in linear algebra and rectangular linear templates [48]. In this context, the rank of a rectangular realvalued template is given by rank of the corresponding realvalued matrix. In analogy to the linear domain, this thesis defines a notion of matrix rank within minimax algebra which leads to a notion of rank of rectangular morphological templates. 20
PAGE 27
21 3.2 Belts and Blogs Consider an algebraic structure (B, o) consisting of a set B together with a binary operation o which assumes the role of addition in B. We call (B, o) a commutative band or semilattice if the following axioms hold for all x,y,z 6 B: Xi : X <> y Â— y o X . X2 ' X o {y 0 z) = {x oy) o z . X3 : X o x = x . Homomorphisms between commutative bands (Bi,o) and (B2, o) are considered to be functions / : Bi Â— > B2, which satisfy the relation f{xoy) Â— f{x) o f{y). Although denoted by the same symbol o, the addition operations in Bi and in B2 are not necessarily the same. Usually, the operations o can be thought of as either the maximum operation or the minimum operation. For these interpretations of o, an additional axiom Xg holds for all a;,?/ 6 B: Xo : X o y Â— x or x o y = y . The axioms Xi, X2, and X^ induce a partial order on B by means of the following definition: x>y<^x = xoy \/x,y^B. Recall that a set B is partially ordered if and only if the following axioms hold for all x,y,z G B: Y\ : X > y and y > x =^ x = y . Y2 : X > y and y > z ^ x > z .
PAGE 28
22 Another partial order on a commutative band B is given by the relation <, where x x. We write x > y or y < x whenever x > y and x ^ y. A commutative band (B, o) satisfying the additional axiom Xq is called linear. In this case, the relation > establishes a linear order on B and satisfies an axiom Yq corresponding to the axiom Xq. Yq : x > y or y > x . The real numbers R as well as the extended real numbers RÂ±oo together with V or A are examples of linear commutative bands. A belt is an algebraic structure (E, o, *) consisting of a set E and two binary operations o and *. The operation o represents the addition in E, and the operation * represents the multiplication in E. The system (E, o, *) satisfies the following axioms (x,y,z e E): Xi,X2,X3 : (E, o) is a commutative band. X4 : X * {y * z) = {x * y) * z . X5 : x * {y o z) = {x * y) o {x * z) . Xq : [y o z) * X = {y * x) 0 {z * x) . Note that the algebraic structure (E, o, *) is known to be a semilatticeordered semigroup. However, since this term seemed so unhandy, CuninghameGreen coined the term belt, evoking a little flavor of the words ring and band. Like in linear algebra, there is a natural way to define homomorphisms between algebraic structures such as belts. Given two belts (Ej, o, *) and (E2, o, *), a belt homomorphism from Ei to E2 is a function / : Ej ^ E2 which satisfies the conditions stated below. The operations of addition and multiplication in Ej are not to be confused
PAGE 29
23 with addition and multiplication in E^. fix<>y)= f{x)0f{y) , f{x*y) = f{x)*f{y) . Va;,2/Â€Ei. If the commutative band (E, o) underlying a belt (E, o, *) is linear, i.e. if this algebraic structure also satisfies axiom Xq, then the belt (E, o, *) is called linear as well. We speak of a commutative belt whenever the operation * is commutative, i.e. whenever the following axiom X7 holds: Xy : X * y = y * X \/ x,y ^ E. A belt satisfying additional axioms Xg and Xg will be called a division belt. Xs Â• 3 0 e E such that 0*x Â— x*0 Â— x VxeE. Xg : y x Â£ E 3 x 6E such that x' * x = x * x' = 0. Axioms Xg and Xg establish the existence of an identity element and of inverse elements with respect to the multiplication *. As usual, the multiplicative inverse of X G E is denoted by x~^. The real numbers together with V and +, denoted by (R, V, +), form a division belt. It is easy to recognize that the following systems belong to the class of linear belts, but not to the class of division belts: i. (R_oo, V, +); ii. (R+00, V, +); iii. (Rl^, V, +); iv. (R+^, V, +).
PAGE 30
24 Lemma 3.2.1. The following properties hold for arbitrary elements x,y,z of a linear division belt (E, o, *). 1. {X o z > y o z X * z > y * z and z * x > z * y 2. x>y=^x*z>y*z and z * x > z * y . 3. z > X and z > y z > xoy . Proof. For a proof of Property 1, we refer to Proposition 22 and Proposition 25 in CuninghameGreen's Minimax Algebra. Suppose that x > y for some y G E. By Property 1, the inequality x * z > y * z holds for all 2 G E. The assumption x * z = y * z leads to an immediate contradiction, if we multiply by z~^ on the right. Hence, x*z>y*zas desired. Similary, we obtain z * x > z * y. Let z > x and z > y for some x,y,z G E. Consider z o{xo y). By the axiom of associativity Xg, zo{xoy) = {z0x)oy. Using the linearity of the belt (E, o, *) and the facts that z > x and z > y, we obtain that z o {x0y) = z, but zo{xoy)^xoy. In other words, z > {x0y).
PAGE 31
25 Blogs are progressively constructed from division belts (E, o, *). First a dual addition o' is introduced for all x,y 6 E. X o' y = [x~^ o y~^) ^ Next we define a set O by adjoining universal bounds Â— oo and +00 to E, i.e. these new elements satisfy (Â—00) 0 X = xo (Â—00) = X (+00) o a; = a; o (+00) = +00 (Â— 00) o' X = a; o' (Â— 00) = Â—00 (+00) o' X Â— xo' (+00) = X Vx G O . Then we extend the domain of the operation * to O x O as follows: (Â—00) * (+00) Â— (+00) * (Â—00) = Â—00 (Â—00) * X = X * (Â—00) = Â—00 V a; G E U {Â—00} (+00) * X = X * (+00) = +00 V ar Â€ E U {+00} . Finally, we introduce a dual multiplication *' which acts exactly like * except that we stipulate (Â—00) *' (+00) = (+00) *' (Â—00) = +00 . Hence, a blog is given by the system (O, o, *,o', *'). Note that the systems (O, o, *) and (O, o', *') form belts. The set E is called the set of finite elements of this blog. If (E, o, *) forms a linear division belt, then we speak of a linear blog. Similarly, (E, o, *) forms a commutative division belt, then we speak of a commutative blog.
PAGE 32
26 By a direct consideration of cases, we can confirm that the following "associative inequalities" hold in an arbitrary blog (O, o, *,o', *'). Lemma 3.2.2. If {O, o, *,o', *') is a blog, then all x,y,z Â€ O satisfy X * (ij *' z) < (x *?/)*' 2 , X *' {y * z) < (x *' J/) * 2 . The principal interpretation is ^RÂ±oo, V, A, +, +'^. The operations + and +' act like the usual addition with the following exceptions. (oo) + (+oo) = (+oo) + (oo) = oo {oo) +' (+oo) = (+oo) +' (oo) = +00 . CuninghameGreen shows that the only finite blog is the threeelement blog representing the set {oo, 0, +00} together with the operations V, A, +, +' restricted to this set. 3.3 Spaces of Vectors and Matrices In analogy to the concept of vector space over a field in linear algebra, CuninghameGreen defines a concept of bandspace, or briefly space, over a belt. Suppose (B, o) is a commutative band and (E, o, *) is a belt. Furthermore, assume that a right multiplication of elements of B by elements of E exists, i.e. there is a binary operation * : B x E Â— > B. Then we say (B, o) is a (right) bandspace over (E, o, *), or briefly B is a space over E, if the following axioms hold for all x,?/ Â€ B and for all A,/z Â€ E: 51 : X * {\ * fi) = (x * \) * fi . 52 : [x 0y) * X = {x * X) 0 [y * \) . 53 : X * (Xo fi) = [x * X) o (x * fi) .
PAGE 33
27 Additionally, if E has an identity element 0 with respect to the multiplication *, we assume for all a; 6 B: : X * 0 = X . Clearly, we can formulate "left" variants of axioms Si to 54 and define thereby a left bandspace over E. We will call B is a 2sided space over E, if B is a right bandspace over E as well as a left bandspace over E, and if the following associative law holds for all a; e B and for all \,n Â€ E: X * {x * fi) Â— {X * x) * fi . For any given belt E and integers m and n, the set ofmxn matrices with entries in E, which is denoted by E"*^", forms a 2sided space over E. In this setting, the addition o in E"*^" as well as the scalar multiplication are performed entry wise. Likewise, the set E" of ntuples with entries in E represents a 2sided space over E. Addition and scalar multiplication are induced by addition and multiplication in E. The algebraic operations o and * can also be used to define a matrix product * : E"*^'' X E^^" + E"*^", where p 6 N. In analogy to the usual matrix product known from linear algebra, multiplying a matrix A e E."^^p by a matrix B Â€ E''^" results in a matrix C Â€ E"*^", whose entries are given by the following formula: Cij = Or=i{aiT *brj) (31) Vi = l,...,m; Vj = . In particular, if E equals the belt (RÂ±oo,V, ), the matrix C representing the matrix
PAGE 34
28 product of A and B is computed as follows: r=l = l,...,m; Vi = l,...,n . In this case, we denote the matrix product by the symbol El . Similarly, the matrix product, which is derived from the operations A and +' in the belt (RÂ±cx>, A, +'), is denoted by the symbol El . A matrix product u * v of a column vector u and a row vector v is called an outer product. We say that a matrix A decomposes into k vector pairs, if A can be written as a k sum of k outer products. For the matrix product El , we have A = V M v'), where u are column vectors and v' are row vectors for / = 1,. . .,k. Example. 8 7 0 9 0 1 I 0(3 8 7) .0 9 0) Theorem (CuninghameGreen) 3.3.1. The following associative inequalities hold for matrices over RÂ±oox^{Ymz) > {x\siY)Enz , XED{Y]SZ) < {XmY)\S\Z yxe R^^" ,v r e R^^Â° yze RÂ±i
PAGE 35
29 Proof. Using distributivity in the belt (RÂ±oo,V, +) and Lemma 3.2.2, we infer for arbitrary i 6 {!,..., m} and j G {1, . . . , p} o n {{X El F) El Z)ij = V A (^'^ +' y^h) + zhj h=l k=l 0 n = V A [i^tk +' vkh) + zhj] h=l k=l 0 n = V A [^"t +' {vkh + zhj)] h=l k=l 0 n The principle of opening implies that V f\ [xik +' [ykh + ^hj)] < h=l k^l 0 V [^ik +' {vkh + Zhj)] for all A; = l,...,n. By applying the principle of h=i closing and distributivity in the belt (RÂ±oo, A, +'), we obtain on no V A [^"t +' (y^h + zhj)] < A V [Â•'"'fc +' (yf'h + zhj)] h=\ k=l k=l h1 A xik +' V {Vkh + Zhj) h=i = (X0(rE]Z)),, . For a proof of the second half of the theorem, we refer to the proof of Theorem 63 in CuninghameGreen's Minimax Algebra. The matrix product defined in Eq. (31) can be viewed as a special case of the generalized matrix product [9]. In conjunction with the matrix sum o and the matrix product *, the set E"^" forms a belt. Direct verification of the relevant axioms shows that E''**" is a commutative band as well as a right bandspace over the belt E"^". Another important class of spaces over a given belt (E, o, *) is the class of function spaces, that is to say spaces in which the commutative band (B, o) is the commutative band (E^, o) of all functions from some given set U to E. Here the addition in E^ and the scalar multiplication on E^ are induced by the addition and the multiplication
PAGE 36
30 in E. Specifically, we obtain the following definitions for all u e U, for all f,g Â€ E^, and for all A G E: {fo9){u) = f{u)og{u) , (/*A)(u) = /(Â«)* A, {X*f){u) = .*f{u). In view of these definitions, it is easy to see that the function space E^ is a 2sided space over E. For example, we can take E to be the real numbers. Although not defined in his book Minimax Algebra, CuninghameGreen suggests that bandspace homomorphisms are functions F from one bandspace B over a given belt E to another bandspace C over E which satisfy the following criteria: F{bi0h2) = F{bi)oF{b2) V6i,62eB F(A * 6) = A * F{b) VAeE, V6GB. The notions isomorphism, endomorphism, and epimorphism have the usual meaning whenever they refer to homomorphisms between commutative bands, belt homomorphisms, or bandspace homomorphisms. If M and N denote the sets {l,...,m} and {1, . . . ,n} for given integers m and n, then the function space e'^^^ can be identified via an isomorphism with the space of m x n matrices over E. 3.4 Conjugacy In analogy to conventional linear operator theory, the theop' of minimax algebra gives rise to a theory of conjugacy. Taking an axiomatic conjugacy approach, CuninghameGreen defines a certain involution x x* subject to certain axioms. In the present section, we discuss axiomatic conjugacy for belts and matrices over belts.
PAGE 37
31 A commutative band B2 is said to be conjugate to a given a commutative band Bi, if Bi is band isomorphic to B2. Suppose now that Ei and E2 are belts. The belt E2 is conjugate to the belt Ei, if there exists a belt isomorphism g from Ei to E2. Note that conjugacy is an equivalence relation, in particular Ei is conjugate to E2 with respect to the isomorphism g~^. Hence, we can say that E2 is the conjugate of Ei. Vice versa, El is the conjugate of E2. We denote the conjugate of a given algebraic structure E by the notation E*. Hence (E*)* = E. Just as in conventional linear operator theory, we shall henceforth employ the notation X* to denote the element of E* corresponding to x e E under the isomorphism from E to E*. Thus (x*)* = X. We call the element x* e E* the conjugate of a: G E. A straightforward verification yields that every division belt is conjugate to itself, or briefly selfconjugate, under the isomorphism given by x* Â— x~^. We say that every blog (O, o, *,o', *') is selfconjugate, because the belt (0,o', *') is the conjugate of the belt (O, o, *) under the following isomorphism: (Â— cx))* = +00 {+00)* = Â—00 ^* whenever x is a finite element of O. Recall that, for any integers m and n, the set of m x n matrices over a belt E forms a commutative band and the set of n x n matrices over a belt E forms a belt. Hence, it is possible to consider the conjugate systems of (E"*^", o) and (E"^", o, *). In fact, Theorem 7.4. of Minimax Algebra shows that the conjugate systems (E"*^")* and (E"^")* are exactly the systems (E*)"^"* and (E*)"^". A matrix A Â€ E"*^" is mapped to its conjugate matrix A* G (e*)"^"* under the
PAGE 38
32 isomorphism of commutative bands from (E'"^", o) to ^(E"*^")*, o^. The matrix A* is obtained from A by transposing and conjugating, i.e. for alH = 1, . . . ,m and all j Â— 1, . . . ,n, the (Â«,i)th entry of A* is equal to the conjugate of the (j, 0th entry of A. Formally, we write: The same relationship as above holds for the belt isomorphism from (E"*^", o, *) to ((E"X")*, o, *). 3.5 The Principles of Opening and Closing The principles of opening and closing provide a set of rules which are useful for the purpose of proving theorems in minimax algebra. Let (W, o, *,o', *') be a blog, and let the symbol V denote the set of finite elements of this blog. Recall the definition of the dual addition o' on the set of finite elements of W: X o' y = [x~^ o y~^) ^ This equation implies that the following equation holds for all xi,X2, . . . ,a;Â„ Â€ V: [oui^sV = WUi^^r' (32) We now discuss some of the basic formal properties of the operators 0 and '. We first address the principle of closing. Suppose that xi,X2, . . . , a;Â„ and j/2, Â• Â• Â• , 2/n are elements of the blog W. Let the following inequality be given: Xi>yi Vi = l,...,n (33) ^ Xi0yi = Xi V i = 1 , . . . , n .
PAGE 39
33 By using the commutativity and associativity of the operations of addition o and dual addition o', we can show that: OU{^^)>0?=l{y^), (o')L(.)>(o')Li(2/.)Since the operation o is idempotent, i.e. x o x = a; for all a; Â€ W, it may be possible to use the following facts in order to simplify the inequalities in Eq. (34): x,= or V Â« = 1, . . . ,n =4> Of^i(a;,) = (O')Li(^') = ^ ' y.= X V i = 1, . . . , n ^ OUiy,) = (^')Li(f') = ^ Â• The manipulation principle used in order to obtain Eq. (34) from Eq. (33) is called the principle of closing. This principle is expressed in Proposition 31 of Minimax Algebra: Theorem 3.5.1. IfW is a blog, then an inequality, governed by a universal quantifier, may be replaced by an inequality having the operator () (respectively <0>'j governing each side of the original inequality, with the index in both cases running over the same finite subset of W as for the universal quantifier. When applying the principle of closing, we add the operator <0> or the operator C)' to both sides of an inequality. CuninghameGreen also introduces the principle of opening, which, informally speaking, relies on the following facts: 1 . An expression governed by an operator 0 is decreased, if the operator is removed. 2. An expression governed by an operator <0>' is increased, if the operator is removed. For example, let us consider the expressions (>^i{xi) and (0')"=i(^Â«)' where xi,X2,. . . ,a;Â„ Â€ W. The associativity and commutativity of the operations o and o'
PAGE 40
34 yields two sets of inequalities: Oi=i{xi)>xi Vi = l,...,n The validity of these inequalities is preserved when applying identical sequences of additions, dual additions, left and right multiplications, and left and right dual multiplications on both sides. These considerations lead to the following theorem: Theorem 3.5.2. Let an algebraic expression /3 be obtained from a wellformed algebraic expression a by deleting an operator 0 which does not form part of the argument of any inversion operation ()~\ Then the inequality a > (3, followed by a universal quantifier whose index runs over the same (finite) set as for the operator 0, is valid in W. Dually, let an algebraic expression ji be obtained from a wellformed algebraic expression a by deleting an operator X, * [OLi(x,)]' V i = 1, . . . , n .
PAGE 41
35 3.6 Linear Independence and Rank In linear algebra, the rank of a matrix is defined as the number of linear independent row vectors or column vectors. In this section, we present CuninghameGreen's definition of matrix rank in minimax algebra, which is based on the concept of strong linear independence. Then we provide our own definition of matrix rank and relate this notion to the one given by CuninghameGreen. As a matter of convenience, we assume that (W, o, *,o', *') is a commutative and linear blog. We denote the latticeordered group of the finite elements of W by the symbol V. Definition. A vector v G W" is said to be linearly dependent on the vectors . . . , v'' 6 if and only if there exist c, G W for alH = 1, . . . , such that v = Oti(c. *v') . We define linear independence as the mere negation of linear dependence. Vectors y^, ... 6 are linearly independent if no one of them is linearly dependent on the others. Example. Let the following elements of RÂ±cx be given: Since v = [2 + vi] V [(Â—2) + V2], the vector v is linearly dependent on and v^. The book Minimax Algebra by CuninghameGreen provides a convenient mechanical procedure, called the ^test, in order to check if given vectors v^, ... , G W" are linearly independent. The 2ltest involves the formation of a matrix 21 6 VV"^". ff o,j
PAGE 42
36 /nxk denotes the (i, j)th entry of 21, and if V e W^^" is the matrix, whose jth column is v' for all j 1, . . . , k, then the matrix 21 is defined as follows: Â—oo else The 2ltest compares each column of V with the corresponding column of V * 21, and comes to the following conclusion: Theorem 3.6.1. For each j = l,...,k, the jth column o/V * 21 is identical with v' if and only ifv^ is linearly dependent on the other columns ofVe W"^'^. The elements of the jth column of% then give suitable coefficients to express the linear dependence. Example. We illustrate this technique with an example from the principal interpretation (rÂ±oo, V, A, +, f'y Let v\ . . . , be the following vectors in The next step consists of forming the matrix V* El V, with v\ ... , v"* representing the column vectors of V: /I 3 oo \ 3 4 5 2 2 5 \3 1 3/ (00 1 1 2^ oo 0 2 3 oo 0 02 Voo 01 0 / 1 3 2 3' 3 4 2 1 Â— oo 5 5 3
PAGE 43
37 When overwriting the main diagonal of this matrix with oo, we obtain the matrix 21 Â€ RÂ±^. We proceed by forming the product V 0 21. /oo 1 1 2 \ _ Â— oo Â— oo Â—2 Â—3 oo 0 oo 2 ' \Â— oo 0 Â—1 Â— oo / /oo 3 2 0\ VElSl = oo 4 2 0 j . \oo 5 3 3/ A comparison of V with V El 21 shows that is linearly dependent on v\ v^, and with coefficients occurring in the second column of 21, i.e.: v= = (1 + v') V (0 + v3) V (0 + v4) . In conventional linear algebra, the concept of linear dependence provides a basis for a theory of rank and dimension. In minimax algebra, however, some unexpected anomalies occur regarding the notion of linear dependence. CuninghameGreen addresses this situation in the following two theorems (Theorems 164 and 165 of Minimax Algebra). Theorem 3.6.2. Suppose W is a blog other than ({oo, 0, oo}, V, +, A, +'). Ifn>2 and k > 1, then there exist k finite ntuples, no one of which is linearly dependent on the others. Theorem 3.6.3. IfW is the threeelement blog ({Â— oo,0,oo}, V, +, A, +'), then there exist (n^ Â— n) finite ntuples, no one of which is linearly dependent on the others. The dimensional anomalies, expressed in the theorems above, lead CuninghameGreen to introduce the idea of strong linear independence.
PAGE 44
38 Definition. Vectors v\ ... , v'' 6 W" are called strongly linearly independent (SLI) if and only if there exists a finite element of W", i.e. a vector v e V" such that v has a unique representation v = 0;^i(a.*v'), a. GW. If y G W"'^^ is the matrix whose columns are v\ . . . , v'' G W", then the strong linear independence of v\ . . . , v'' G W" is equivalent to the following condition: There exists a vector v G V" such that the equation V * x = v is uniquely soluble. Another equivalent condition is given by the theorem below, whose proof takes advantage of the following lemma: Lemma 3.6.4. Suppose v\ . . . , v'' G V". The vector v = Of^j (Ov') ^ * v' G V" is bounded from above by the ndimensional vector, whose entries are all equal to ly, ic. vj <1\/ V j = 1, . . . ,n . Proof. Let us consider the algebraic expression given by (OJLi^^/) ^ * ^]=i'^j Opening with respect to the index j yields the following inequality, which is valid in V Oil {otxv\)~'*v]
PAGE 45
39 Proof. "SLI Eq. (35)": If given vectors \\ ... 6 V" are SLI, then v = =iv' < <>,=iv' for all p = I,. . . ,k imply that Eq. (35) holds. "Eq. (35) =^ SLI ": Given vectors v\ . . . , e V", we assume that Eq. (35) holds, i.e. 0Liv' < 0Liv' for all p = 1, . . . , A;. We set a, = (o^^^v()~'= (Ov')"'We will show that the vector v = 0 Lj (a^ * v') G V" does not have another representation. Hence, if v = 0 Lj (c, * v*), it suffices to show that c, = a, for alH = 1, . . . , k. The following argumentation shows that c, < a, for all ?' = 1, . . . , k. The statement V = OiLi {ci * v') is equivalent to Vj = Of^i (ci * j for all j = 1, . . . , n. For fixed, but arbitrary j 6 {!,...,Â«}, consider the algebraic expression Of=i (ci * v*^^ . An application of the principle of opening with respect to the index i yields: c, * v'j < <;^^=i{ci * v)) yi = l,...,k. By assumption, <)f=i (ci * u]^ = vj for all j = 1, ... , n, and by the previous lemma the entries vj of v are bounded from above by lyTherefore, Ci*Vj<\\/ y i = I, . . . ,k; V j = 1, . . . , n . A right multiplication by (v'j^ yields: Ci<{v)y" V^ = 1,...,A;; Vi = l,...,n. We are able to finish the proof of c, < a, for alH = l,...,k by closing with respect
PAGE 46
40 to the index j and by using Eq. (32): Assume that there exists apG {!,..., n] such that Cp < ap. Then Cp*v^ < ap* for all j = 1, . . . ,n. By the principle of opening, ai*v*j < 0,=i(ai *t;]) = vj Vi = 1,...,A;; Vi = . Thus, Cp * < vj for all j = 1, . . . ,n, which implies that Cp * < v. Combining the inequalities Of=i(c, *v') < v and Cp * < v, where j = l,...,n, leads to a contradiction to the assumption v = Of=i(ct * v'): Of=i {ci * v') = Of=i(c, * V*) O (Cp * vP) < V . Example. Let the following elements of be given: Note that Â• the corresponding normed vectors v, , e = 1,2,3 are given by 2^ vi 0 V2 = [si {2 Â• The vectors vi, V2, vs are not SLI since 3 VI V V2 = I 0 j = Y V, .
PAGE 47
41 The following theorem shows that the concept of strong linear independence is more suitable for defining notions of rank and dimension in minimax algebra than the concept of linear independence. Theorem 3.6.6. For any linear blog W there are k vectors v\ ... , G W" which are SLI if and only if k G {1, . . . , n). Definition. We define the SLIrank of a finite matrix A G V"*^", denoted by ranksLi(^), as the maximal number of row vectors which are SLI. Theorem 3.6.7. The SLIrank o/ A G V"*^" equals the maximal number of column vectors of A which are SLI. Note that by Theorem 3.6.6 the SLIrank of a finite matrix A G V'"'*" is less than or equal to min {m, n}. In linear algebra, the concept of matrix rank naturally induces a concept of rank of rectangular templates. This notion has found many applications in the area of template decomposition due to the fact that it represents the minimal number of those products of row and column templates which add up to the original template. These considerations lead us to introduce a concept of matrix rank within minimax algebra which has similar consequences with respect to the purpose of template decomposition. Definition. We define the rank of a finite matrix A G W"*^" as the minimal number of vector pairs into which it can be decomposed, i.e. the minimal number r for which there exist column vectors u\ . . . ju"^ G W"*^^ and row vectors w^, ... , w"^ G W^^", such that A has a representation of the following form:
PAGE 48
42 Theorem 3.6.8. Let W be a linear blog. The rank of a finite matrix A Â€ W"^" is bounded from below by the SLIrank of A and bounded from above by the maximal number t of linearly independent row vectors of A as well as by the maximal number of linearly independent column vectors of A. Proof. "ranksLi(^) < rank(yl) "; Recall that the rank of A is defined as the maximum number k of SLI row vectors. Let v\ . . . , 6 W^^" be A; SLI row vectors of A. By Theorem 3.6.5 Suppose that rank(yl) = r. Then there exist row vectors d\...,d'" G W^^" such that each row vector of A can be expressed as a linear combination of d\ . . . , d*". In particular, each row vector v*, where i = l,...,k, has such a representation. Hence, Let c/ = <>/=i (c) and let v = <0>f_jv'. Using the axioms of associativity, commutativity, and distributivity, we show that v = 0}"_i(c/ * d'): OLiv'
PAGE 49
43 Given any / Â€ {!,... ,r}, there exists and index z/ G {1, . . . , A;} such that q = c' due to the linearity of W. This fact yields q * < v" in light of the equality ^cj' * d' j = v". Closing this inequality with respect to the index / yields 0[=i(c/ * d') < O/^jV". Thus, we are able to conclude the proof of the inequality rank(A) > ranksLi(^) as follows: V = oui (ci * d') < oLiv" < oly = V =^ [=lV" = V =^ \{h,...,ir} \ > k =i~ r>k . "rank(y4) < t ": It suffices to show that there exist column vectors u\ . . . ,u*' e W"^^ and row vectors w\ ... , e VV^^", such that A has a representation of the following form: A=0U(u'e1v') . Let a(z) G W^^" denote the ith row vector of A for all i = l,...,m. WLOG the first t row vectors of A are linearly independent, so each row vector a(z) can be expressed as follows in terms of a(l), . . . ,a(^): a(0 = (c) * a(i)) Vi = l,...,m, where cj G R are some appropriate scalars for all i = 1, . . . ,m and for alH = 1, . . . , /. If c' G R"*^! denotes the column vector V/=l,...,f,
PAGE 50
44 we have A = OLi(c'*a(/)), which concludes the proof of the theorem.
PAGE 51
CHAPTER 4 RANKBASED DECOMPOSITIONS OF MORPHOLOGICAL TEMPLATES 4.1 Introduction By way of a canonical bijection, matrices correspond to rectangular templates. Hence, the notion of matrix rank induces a notion of rank of rectangular templates both in the linear domain and in the nonlinear domain. Moreover, algorithms for the decomposition of matrices translate into algorithms for the decomposition of templates. A rankbased decomposition of a rectangular template is considered to be a weak decomposition of a rectangular template according to its rank, i.e. a decomposition into a minimal number of sums of outer products of row and column templates. Rankbased decompositions have the following benefits: The computational complexity of a convolution of an image with a rectangular template reduces from a problem which is quadratic in complexity to a linear problem. The rankbased decomposition of a template can be used as a first step for the decomposition into a form which matches the neighborhood configuration of a particular parallel architecture. The strong decomposition of a rank 1 template is an easy task both in the linear and in the nonlinear domain [38]. O'Leary [44] showed that any linear template of rank 1 can be factored exactly into a product of 3 x 3 linear templates. Templates of higher rank are usually not as efficiently decomposable. However, the LU factorization yields a proven method for determining a rankbased decomposition of a linear template of arbitrary, unknown rank [28, 48]. 45
PAGE 52
46 Thus far, rankbased decompositions of morphological templates are restricted to decompositions of rank 1 templates. In this chapter, we provide a polynomial time algorithm for the rankbased decomposition of morphological templates of rank 2. Furthermore, we present a heuristic algorithm for the rankbased decomposition of morphological templates of arbitrary rank. 4.2 Correspondence between Matrices and Rectangular Templates Note that there is a natural bijection (f> from the space of all m x n matrices over R_oo into the space of all rectangular m x n templates in (R?oo) Â• Let y = G X be arbitrary and x = {xi,X2) G X be such that . fmm Â— l1 .fnn Â— ll xi = X Â— mm < Â— , Â— Â— >, X2 Â— y Â— mm < Â— , , 2 ' 2 J ' ' " 1,2' 2 J The image of a matrix A G R'^x" under (f) is defined to be the template t e (R^oo) which satisfies ty(a;i + i I, X2 + j I) aij V t = 1, . . . , m, V j = 1, . . . n , ty(j/lÂ» 2/2) = 00 V t/i,y2 0 [xi,xi + m 1] X [x2,X2 n 1] . The natural correspondence between rectangular templates in t G (R00) matrices over R_oo allows us to use a minimax algebra approach in order to study the weak decomposability of rectangular templates into separable templates. Example. Let A G R^^^ be the matrix and u, v the vectors given below. 2 5 10\ /2' A = \ i 7 12 , u= 0 , v = (7,4,12) 3 0 5 / V 7
PAGE 53
47 The function {A) for some realvalued matrix A, then we define the rank of the template t G (R?oo) to be the rank of the matrix A. Our interest in the rank of a morphological template is motivated by the problem of morphological template decomposition since the rank of a morphological template t G (R^oo) represents the minimal number of separable templates whose maximum is t or, equivalently, the minimal number r of column templates r' G (R'^oo) and row templates s* G (R'^oo) such that t = V (r'0s'). Lemma 4.2.1. If A e R"*""" andu^v < A, where u G R"*""' and v G R'""", then < A ( ^ij ~ ^! ) f^'^ ^ Â— 1 ) Â• Â• Â• , Proof. Each entry aij of A is greater than or equal to ui + vj for all / = n implying that ui < f\ {aij Â— vj) for all Â« = 1, ... , m. k Theorem 4.2.2. If a matrix A G R"*^" has a representation A Â— \/ (u' M v') in terms /=i of column vectors u' G R"*^* and row vectors v' G R^^", where I = l,...,k, then A can be expressed in the following form: k A^y (w'eIv'), /=1 where w' G R""^^ is given by n j=i
PAGE 54
48 Proof. Let W = \/ (w' 0 v') . The proof is accomplished in two parts. 1. "W < A": = V(w'^v'),^. k V A . \s=l k < V [""'J ~ /=i Â— a, J , i = l,...,m , = l,...,n IE 2. 'W > A": Since A G R"*^" can be represented as V (u'EI v'), we observe that 1=1 by the previous lemma u< /\ ( a,j j for all i, j, and /. The latter inequality has j=i ^ the following consequences, which conclude the proof of the theorem: k /=1 ls=l m K s=l J 1=1 Remark. Theorem 4.2.2 implies that, for any matrix A G R'nx" of separable rank it, it suffices to know the row vectors e R^^", I = 1 . . ,k, which permit a weak decomposition of A into A; separable matrices in order to determine a representation of
PAGE 55
49 A in the form k w' G R V /= l,...,k . Like most of the theorems established in this dissertation, Theorem 4.2.2 has an obvious dual in terms of column vectors which we choose to omit. ^ 4.3 Invariance Properties of the Rank of Realvalued Matrices In this section, we derive some invariance properties of the matrix rank which translate directly into results about the rank of rectangular templates. These theorems greatly simplify the proofs of the decomposition results which we will present in the following sections. Definition. Let A 6 R"*^" and /? be a permutation of {1, . . . , n}. The associated column permuted matrix Pc{A) of A with respect to p is defined as follows: Similarly, if p is a permutation of {1,. . . ,m}, then we define the associated row permuted matrix Pr{A.) of A with respect to /) as follows: pM) = /Â«lp(l) Â«lp(2) a2p(l) Â«2p(2) Â«lp(n) \ Â«2p(n) m, .p(l) Â«mp(2) Â• Â• Â• Omp(n) ' Pr{A) = /ap(l)l ^1)2 Â«p(2)l Â«p(2)2 Â«p(l)n \ ap(2)n p(m)l "p(m)2 p(m)n I Theorem 4.3.1. The rank of a matrix A G R"*^" is invariant under column and row permutations.
PAGE 56
50 Proof. Suppose the matrix A e R'wx" has rank r. Hence, there exist column vectors u' G R'^xi and row vectors v' e R^^", such that A can be written in the form Let /) be a permutation of the set {!,..., n}. An application of the column permutation pc to the matrix A yields the following identities for all i = 1, . . . ,m and for all j = 1, . . . , n: r r ft /=i /=i *^ \l=i Hence, Pc{A) = y u^M Pc{^') , which implies that the rank of Pc{A) is at most r. 1=1 Now it suffices to show that pc{A) can not be weakly decomposed in the form it . , pc{A) Â— V 0 , where G Rmxi, v' Â€ Rixn, and A; < r. Suppose, we are given /=i k . a representation Pc(^) = \/ u^M^^ Applying the inverse of the column permutation Pc to both sides, we obtain the following equation: A=\/u^mp;'{y^) . /=i This equation concludes the proof of the invariance of the matrix rank under column permutations, because k must be at least as big as r = rank(y4). The invariance of the matrix rank under row permutations can be proven in a similar fashion. Recall from elementary linear algebra that the set of realvalued m x n matrices together with the operations of addition and scalar multiplication forms a vector space over R. The operation of addition in R"*^" is performed pointwise and denoted as usual by the symbol +. This addition is not to be confused with the operation V, the addition
PAGE 57
51 in the band space R"*>^". Likewise, we differentiate between the scalar multiplications of the vectorspace R'"^^" and the scalar multiplications of the band space R"***". Theorem 4.3.2. The following operations of the vectorspace R""^" preserve the rank of a matrix A Â€ R""^" in minimax algebra: a. additions of matrices of rank 1; b. scalar multiplications. Proof. Let r be the rank of the matrix /I e R"*^". By definition of the matrix rank, the r integer r is the minimal number such that A = V (u' El v') for some column vectors /=i G R"*^^ and row vectors v' G R^^". Suppose, that we are given a separable matrix B G R"*^" and a scalar c G R. The separability of B implies that B can be written as an outer product of a column vector u G Rmxi and row vectors v G RixnThe validity of the theorem follows from the following facts, which are easily derived: r A + 5 = V [(u' + u) 0 (v' + v) r c Â• A = Y Â• u' 0 c Â• v') . 1=1 4.4 Separable Matrices and Templates At this point, we are finally ready to tackle the problem of determining rankbased decompositions of matrices in minimax algebra. The reader should bear in mind the consequences for the corresponding rectangular templates. For any matrix A G R'^^^", we use the notation a(Â«), where i = 1, ... ,m, to denote the ith row vector of A and we use the notation a[j], where j = 1, . . . ,n, to denote the jth column vector of A.
PAGE 58
52 Theorem 4.4.1. (Li and Ritter [38]) Let A Â€ R"'''"i>ea separable matrix and {a(i) : I = 1, . . . , m} be the collection of row vectors of A. For each arbitrary row vector a(io), where I < h < m, there exist scalars Â£ R, i = 1, . . . ,m, such that the following equations are satisfied: A,+ a(io) = a(i), V i = 1, . . . , m . In other words, given an arbitrary index I
PAGE 59
53 STEP201) + 1, j = 2,...,n: Compare ci with aij Â— a\j. If c, 7^ a,j Â— a\j, the matrix is not separable and the algorithm stops. If step 2(n Â— 1) + 1 has been successfully completed, then A is separable and A is given by c21a(l). Note that this algorithm only involves (m Â— l)n subtractions and (m Â— l)(n Â— 1) comparisons. Hence, the number of operations adds up to 2(m Â— l)n Â— m + 1, which implies that the algorithm has order 0{2mn). As mentioned earlier, the given image processing hardware often calls for the decomposition of a given template into 3x3 templates. In the case of a separable square template, this problem reduces to the problem of decomposing a column template into 3x1 templates as well as decomposing a row template into 1x3 templates. Theorem 4.4.3. Let t be a square morphological template of rank 1, given by t = rMs, where r is a column template and s is a row template. The template t is decomposable into 3x3 templates if and only ifr is decomposable into 3x1 templates and s is decomposable into 1x3 templates. Proof. "<= ": Let r\ . . . , be 3 x 1 templates such that r = r' 13 . . . El r*^, and s\ . . . , s^" be 1 X 3 templates such that s = s' El . . El s^. The template t decomposes strongly into the 3 x 3 templates r* Ms\ i Â— I, . . . ,k, for the following reason: t = r E s
PAGE 60
54 "=^": Let t = rEs be strongly decomposable into A: 3 x 3 templates. We show by induction on the number k that r is decomposable into 3x1 templates and s is decomposable into 1x3 templates. If A; = 1 then t is a 3 x 3 template and therefore, r is a 3 x 1 template and s is a 1 x 3 template which leaves nothing to show. For the remainder of the proof, we assume that if a square template t = r 0 s, where r is a vertical and s is a horizontal strip template, decomposes into less than k 3x3 templates, then r is strongly decomposable into 3x1 templates and s is strongly decomposable into 1x3 templates. Let a template t of rank 1 be given which can be represented in terms of 3 x 3 templates a', i = l,...,k. t = a' M . . . El . Let us define the template a as a^ El ... El a*^~\ and the template b as a*^. Furthermore, we introduce the following convenient notations for any rectangular pxq template. Recall that for any matrix C Â€ R''^' the symbols c{i) and c[j] denote the ?th row vector of C, the jth column vector of C respectively. The indices i range from 1 to p and the indices j range from \ to q. Let us introduce the following convenient notations for arbitrary rectangular templates c. c\j] = ir\c)[j]) , cbi. = (r'(cL/i)). . Obviously, t = aElb. Also note that the row template t(l) and the column template t[l] can be computed as follows: t(l) = a(l)Elb(l) , t[l] = a[l]Elb[l] .
PAGE 61
55 Since t is of rank 1, all row vectors (/!Â»~^(t(z)) are linearly dependent on <^~^(t(l)). Hence, there exist constants c,, i = l,...,m such that t{i) = c, + t(l). In fact, c, = t{i)i Â—t{l)i for alH = 1, . . . ,m, which implies that t has another representation besides t = rEl s, namely t = (t[l] El t(l). We conclude that r,+s = ^[l],+t(l) for all i = 1, . . . ,m, in particular for i = 1. Thus, s and t(l) only differ by a constant c, and for similar reasons r and t[l] differ by a constant d. s = t(l) + c , r = t[l] + d. We are now able to finish the proof by substituting a(l) El b(l) into t(l) and a[l] El b[l] into t[l]. s = t(l) + c = [a(l) 0 b(l)] + c = a(l) El [b(l) + c] , r = t[l] + d= (a[l] El b[l]) + d = a[l] El (b[l] + d) . Example. Let A be the realvalued 5x5 matrix given below. A = /3 7 4 2 1 \ 11 15 12 10 9 5 9 6 4 3 6 10 7 5 4 \12 16 13 11 10/ = uEl V, where u = 8 2 3 \9/ , V = (3 7 4 2 1 ) The template t = (j){A) is not decomposable into two 3x3 templates since the template r =
PAGE 62
56 4.5 Matrices and Templates having Rank 2 Lemma 4.5.1. If a matrix A G R"Â»x" has rank 2, then there exits a transform T Â— consisting of only row permutations, column permutations and additions of separable matrices Â— as well as vectors u 6 R"'><^ and v G R^^" such that T{A) can be written in the following form: r(A) = o^xÂ» v(u0v) , where and Omxn denotes the m x n zero matrix. Proof. The matrix A G R"*'*" can be represented by means of some column vectors u\ G R"*^^ and some row vectors G R^^": A = (u' 0 v') V (u^ia v"*) . li B = (U*) El (v*) then we compute the matrix A\B Â£ R^^" as follows: A + B = [{n' M v') V (u= M v^*)] + B = [(u^Elv^) + 5]V[(u=0v=) + 5] = 0Â„xÂ„V[Ku^)El(v^v^)]. Let a denote the column vector and let b denote the column vector v\ The entries of this vectors can be ordered by a permutation p : {1, . . . , m} {1, . . . , m}
PAGE 63
57 such that ap(l) < 0,p{2) < < Â«p(m) Â• Applying the row permutation pr to A + B yields the following result: Pr{A + B) = /3r(0mxn V[aElb]) = />r(0mX7i) V/9r(aElb) = Omxn V [pr{a) M b] . Hence, if we introduce the notations u for /9r(a), v for b and T for the addition of B followed by pr, we are able to conclude the proof of Lemma 1: T{A) = 0 V (uElv) , Ul < U2 < . . . < Um Lemma 4.5.2. Let T one of the rank preserving transform of Theorems 4.3.1 and 4.3.2 and A E R"Â»x", jfig transform T maps row vectors of A to row vectors ofT{A), and column vectors of A to column vectors ofT{A). Proof. Pick two arbitrary entries a,;t and an on the ith row of A and observe where they are mapped under row permutations, column permutations, additions of row vectors, additions of column vectors, and scalar multiplications. Theorem 4.5.3. A matrix A G R"*^" has rank 2 if and only there are two row vectors of A on which all other row vectors depend (linearly).
PAGE 64
58 Proof. Since A e R'"><" has separable rank 2 there are vectors u^, 6 R'^^S and vectors v\ G Ri><" such that A can be expressed as follows: A = (u^EIv^) V(u^Elv=') . By Lemma 4.5.1, there is a transform T and vectors u G R"*^^ and v 6 R*^" such that r(A) = o^xnV(uEiv), Ul < U2 < . . . < Um Let B denote T{A). Note that for every index i G {l,...,m} the ith row vector b(^) is given by h{i) = Oixn V (u, +v). The ascending order of the scalars Â«j, i = l,...,m, induces an ascending order of the vectors h{i), i = l,...,m. Using the fact that ui < U2 < . . . < Um, we can also show that h{i) is linearly dependent on b(l) and b(m): b(0 = 0,xÂ„ V(ui + v) = Oxxn V[(U1 +V) V(U.+V)] = Oixn V [(ui + V) V (ui + Wm + V)] = [Oxxn V {Ui + V)] V [{ui Um) + (Wm + v)] = b(l) V [{ui Um) + b(m)] . These relationships between the vectors h{i) and the vectors b(l) and b(m) imply that B can be written as follows: B = [Omxx M b(l)] V [(u Um) M b(m)]
PAGE 65
59 An application of the inverse transform yields the following: A = T\B) = [T'{Onx.)MT\h{l))] V [T\uUm)MT\h{m))] . Recall that, by Lemma 4.5.2, the vectors r~^(b(l)) and T~^{h{m)) are row vectors of A. Hence, there are indices A;, / G {l,...,m} and vectors c% 6 R"*^^, such that A = [c^l2la(A;)] V [c^Ela(/)]. In other words, every row vector a(i), where i G {l,...,m}, is linearly dependent on the same pair of row vectors a(/j) and a(/). Therefore, we have completed the proof of the sufficiency part of the theorem. Obviously, the other direction of the theorem holds. Remark. By Theorem 4.5.3, a matrix A 6 R"*^" has rank 2 if and only if there exist two row vectors of ^4 Â— a(o), a(p) where o, p E {1, . . . , m} Â— which allow for a weak decomposition of A. In this case, an application of Theorem 4.2.2 yields the following representation of A: A=[um a(o)] V [v 0 a(/))] , where u, V Â€ R""^^ , m Ui = /\ {aij Ooj) (i = 1, . . . , m) , j=l m vi = /\ [aij Upj) (i = 1, . . . , m) . Hence, in order to test an arbitrary matrix A G R"*^" for weak decomposability into two vector pairs, it is enough to compare A with {h[i] M a(i)) V (b[j] E a(j)) for all
PAGE 66
60 indices i, j Â€ {1, . . . , m}, where the matrix 5 Â€ R"'^'" is computed as follows: n bij = /\ {ais ajs) s=l Vt = l,...,m, Vj = l,...,m. Algorithm 2. Assume a matrix 6 ^mxn ^geds to be decomposed into a weak El product of 2 vector pairs if such a decomposition is possible. Considering the preceding remarks, we are able to give a polynomial time algorithm for solving this problem. For each step, we include the number of operations involved in square brackets. STEP 1: Compute B G R"'><"'. [m^n subtractions; m^(n Â— 1) comparisons] STEP 2: Form C* = h[i] M a(i) for all i = 1, . . . , m. [m(mn) additions] STEP 3: Form C' V C' for alH, = 1, . . . , m and compare the result with A.]fCÂ°yCP=A for certain o, p= 1 , . . . , m then the algorithm stops yielding the following result: ^ = (b[o] El a(o)) V (b[p] EI a(p)) . If C" V < for all i, j = 1, . . . , m, then /I does not have a weak decomposition into two vector pairs. [At most 2 ^ (mn) comparisons] This algorithm involves at most a total number of m^(3n Â— 1 + (m Â— l)n) operations which amounts to order 0{m^n).
PAGE 67
61 Example. Let us apply Algorithm 2 to the following matrix A e R^^s A = /9 5 5 6 5\ 7 4 4 3 4 12 7 7 9 7 Vll 3 5 8 3/ We compute matrices B Â€ R^^^ and C G R4X5 for all Â« = 1, . . . ,4. B = C4 = b[4] M a(4) /o 1 3 2\ 3 0 6 5 2 3 0 1 ^2 1 4 0 ) /9 5 5 6 5\ a(l) = 6 2 2 3 4 11 7 7 8 7 1 3 3 4 3/ /8 5 5 4 5\ a(2) = 7 4 4 3 4 10 7 7 6 7 1 U 3 3 2 3^ /9 4 4 6 4\ a(3) = 6 1 1 2 1 12 7 7 9 7 1 3 3 5 3^ 1 3 6 1 \ (4) = 6 2 0 3 2 12 4 6 9 4 3 5 8 3 / Comparing the matrices C VC' with A for all ?, j = 1, . . . , m reveals that A = VC^. Thus, A = (b[2]Ela(2)) V(b[4]Ela(4)) Example. Let A e R9^9 be the following matrix, which constitutes the maximum of a
PAGE 68
62 matrix in pyramid form and a matrix in paraboloid form. /18 14 10 6 2 6 10 14 18\ 14 9 4 1 2 1 4 9 14 10 4 1 4 6 4 1 4 10 6 1 4 7 10 7 4 1 6 2 2 6 10 14 10 6 2 2 6 1 1 4 6 7 4 1 6 10 4 1 4 6 4 1 4 10 14 9 4 1 2 1 4 9 14 V18 14 10 6 2 6 10 14 18/ Since matrices in paraboloid and in pyramid form are separable, algorithm 2 should yield a weak decomposition of A in the form [u El a(o)] V [v El a(p)] for some vectors u, V e and some indices o, p Â€ {1, . . . ,9}. Indeed, if 5 Â£ R^^^ denotes the matrix computed by algorithm 2, A = [b[l]Ela(l)]V[b[3]Ela(3)], where b[l] = /11\ 4 5 8 0 12 3 16 , b[3] = 10 12 3 8 0 4 5 \o/ \nJ a(l) = (18 14 10 6 2 6 10 14 18), a(3) = (10 4 1 4 6 4 1 4 9). This weak decomposition of A can be used to further decompose the square templates {h[3] M a(3)) into 3x3 templates. By Theorem 4.5.3, the templates (b[l] El a(l)) and
PAGE 69
63 column templates (b[l]) and (a(l)) and (f){a{3)) are decomposable into 1x3 templates. It is fairly easy to choose 3 X 1 templates r* G (R?oo) and 1 x 3 templates s' e (R?oo) , i = 1 , Â• Â• Â• , 4, such that (l>{h[l]) = ((r^l21r^)Mr3)Elr4 and (b[l] El a(l)) in the form below. (i>{h[l] M a(l)) = [((r' El r^) M r^) El r^] El [((s' El s^) El s^) E s^] . By rearranging the templates r' and s' for z = 1, . . . , 4, we can achieve a decomposition of {h[l] E a(l)).
PAGE 70
64 In a similar fashion, we are able to decompose the template (j){h[S] M a(3)) into four 3x3 templates. Unfortunately, rankbased decompositions of matrices having rank A: > 3 can not be obtained by employing a similar strategy as in Algorithm 2, since a theorem similar to Theorem 4.5.3 does not hold for this class of matrices. This fact is expressed in the following theorem. Theorem 4.6.1. For every natural number k > 3 there are realvalued matrices which are weakly El decomposable into a product of k vector pairs, but not all of whose row vectors are linearly dependent on a single ktuple of their row vectors. Proof. Let " = ^'^y given G N. In order to prove this theorem we show that a matrix A G R"^^" satisfies the statement of the theorem if it has the following properties: Each row vector a(?) of A has two locations ii, 12 G {1, . . . , A:} having value 1. All other entries of a(i) are 0. li i ^ j, then a(z) 7^ a{j). Note that a(i) = v'^ V v*^ for all Â« = 1, . . . , n, where v' G is defined as follows for all / = 1, . . . , fc: Equivalently, the matrix A has representation A = \/ [u^M v') , where the column 4.6 Matrices and Templates of Arbitrary Rank k vectors u' G R"^^ (/ = 1, . . . , fc) are determined by
PAGE 71
65 Therefore, A satisfies the first requirement of the theorem. Let (a\...,a*^) denote an arbitrary A;tuple of row vectors of A. (Note that the index i in a' does not necessarily correspond to the row index.) Assuming this /jtuple allows for a weak decomposition of A into k vector pairs, there exist column vectors k W e R"^^ (/ = l,...,k) such that A can be written as V (^^ El a'). In other words, k a(O = V(^' + ^0 = Â• (41) Now let us pick an index i Â€ {!,...,Â«} such that a(i) 7^ a^ for all / = 1, . . . , For notational convenience, let us use the abbreviations a for a(i) and 6' instead of b\ for all i = 1, . . . , n and for all / = 1, . . . , A;. Since a differs from every vector a' among a\ . . . ,a*^, there has to be an index /' for each / = l,...,k such that ai> = 0 and a\, = 1. Since 6' + a' < a, the scalar 6' cannot exceed 1. These considerations lead to the following sequence of inequalities contradicting the assumption presented in Eq. (41): k k /=i /=i Hence, the row vector a of A is linearly independent from the A;tuple (a^, . . . ,a*). Corollary 4.6.2. Let A; > 3. There exist integers m and n, and matrices A, B e R"*^" such that rank(J9) = k and such that for any metric d on R"*><" rf^^,\/u'Ela(oj >d{A,B) V/C{l,...,n} : \I\ = k; V u' G R"* .
PAGE 72
66 Proof. The previous theorem implies that for every k > 3, there exists a realvalued matrix A satisfying We will show in chapter 6 that the general problem of determining the weak decomposition of matrices is NPcomplete. However, the NPcompleteness of a certain problem does not preclude the existence of an efficient algorithm for solving arbitrary instances of this problem. For matrices having relatively small rank compared to their size, we suggest the following heuristic algorithm. This heuristic is based on the following observations: Remark. Suppose that A G R"Â»^", whose rank is an unknown integer r. Hence, A can be represented as the maximum of matrices A' e R"*><", where / = l,...,r and a' = x' El y' for some realvalued column vectors x' of length m and some realvalued row vectors of length n. According to Theorem 4.4.1, the separability of the matrix a' induces the equality of the differences a\j^ a\j^ for all z = 1, . . . , m and for arbitrary, but fixed / G {1, . . . ,r} and ji, j2 Â€ {1, . . . ,n}. On the other hand, if a,jj a.jj for i G / C {!,..., m}, it seems reasonable to assume that Oij^ Â— a\j^ = a;' j/j^ and a,j2 = a\j^ = x'+ y'j^ for most i G / provided that "r = rank(/l) is relatively small and that {a,j : i = 1, . . . , m; j = 1, . . . , n} is relatively large compared to r". This assumption and the following two lemmas play an important role in Algorithm 3 which intends to determine r Â— ra.nk{A) as well as column vectors u' G RÂ™*^^ and row vectors V/C{l,...,n} : 7 = A:; V u' G R m
PAGE 73
v' Â€ Ri^'^ such that 67 /=1 Lemma 4.6.3. Let A e R'"''" and x,u Â€ R""""^ and y,v e R^^". Suppose that there exists an index j' G {l,...,n} such that + yji = aiy = + Vji for all i G / C {1, . . . , m}. i/"/or all j G J ^ {1, . . . , n} there exists an index ij G / such that + Vj = ci^j > + yj, then Ui + Vj > Xi + yj V I G /, V J G J . Proof. Ui + Vj {ui W.j) + (Wtj + Vj) = [{ui + Vj>) (u, . + Vj')] + {Ui^ + Vj) = [aij' aijj'] + {ui, + Vj) > [{xi + yj') {xi, + yj')] + {xi^ + yj) = [xi X, .] + (a;, . + yj) = Xi + yj yiei^Wj gJ . Lemma 4.6.4. Suppose u,x G R'"^\ v,y G RL^, and jij2 e {1, . . . ,m}. Furthermore, let a subset / o/{l, . . . , m} and a subset J of {1, . . . ,n} be given such that vj and Xj are finite for all j G J. If the following conditions are satisfied then ui + Vj > xi + yj
PAGE 74
68 for all i = 1 , . . . , m and for all j G J. ui + = Xi + = G.jj y i e I , Ui + Vj^ = Xi + yj^ = Cij^ y i E I , Â«t' + ^ji = or Ui' + = aij2 V Â«' Â€ {1, . . . , m} \ / , Ui + vj > Xi + yj MiEl^yjEJ. Proof. It suffices to show that for arbitrary i' Â€ {1, . . . , m} \ / and for arbitrary j G J the relation Uii + Vj > x,/ + yj holds. Suppose / denotes ji if + Vj^ = a,jj and j2 otherwise. We obtain the following estimate for the difference u, Â— u,', where i is an arbitrary element of /: Uii Ui = {ui' + Vji) {Ui + Vj') = ai'j' Â— Giji > {xi' + yj') {xi + yj') Using this estimate, the inequality u,' + Vj > xi' + y^ now follows from the inequality Ui + Vj > Xi + yj. Remark. Assume that A = \/ (x'My'), where x' and y' are unknown realvalued /=i vectors. Algorithm 3 constructs vectors u' and v' for / = 1, . . . , r such that u' 0 < A. k k The maximum V (u' El v') reaches A, i.e. V (u' El v') = A, if there exist sets and /=i /=i
PAGE 75
69 for all / = 1 , . . . , such that I xj' C {(?,;) : xi + yj = aij} I X j' C {{i,j) : Xi + yj < ui + vj} Lemmas 4.6.3 and 4.6.4 already provide some insight into these matters. After presenting Algorithm 3 and applying this algorithm to a specific example, we establish a theorem which should lead to an even better formal understanding of Algorithm 3 by combining Lemmas 4.6.3 and 4.6.4. Algorithm 3. Let A e R"Â»x". This algorithm determines a parameter k > rank(A) and progressively constructs column vectors u' G R"*^^ and row vectors G R^^'^ such that Furthermore, this algorithm generates a matrix M Â€ R'"'*" which has the following property: Initialize the matrix M 6 R"*^" to be the zero matrix Omxn and let / = 0. As long as the matrix M has a zero entry, we increment the parameter / by 1, and perform the following steps. Otherwise k = I and the algorithm stops. k = 1, . . . , m; V; = 1, . . . , n; V / = 1, . . . , A;
PAGE 76
70 STEP 1: (Find the most frequent tuple difference) Determine a strictly increasing sequence of maximal length (ii, 12, Â• Â• Â• , ip) such that there exist two different column indices ji and j2 having the following property: and ^itji = 0 or rui^j^ = 0 Vi = l,...,p. For any t = satisfying m,,j, = mi^j^ = 0, we set m,,_;i = m,,j2 = /. If either one of the entries mi,j^ or mi^j^ differs from 0 for alH = 1, . . . ,p, then set mjjj, = ruiij^ = / for alH = 1 , . . . , p. STEP 2 (Extend the markings of found tuples horizontally): For all i Â€ . . . , ip} and all j = 1, . . . , n, we set m,j = / if a. m,jj /. m b. Uij aij^ = /\ {a^jazj^). 2=1 STEP 3 (Create preliminary row vectors v' G R!",^'): Initialize the vector v' e H^^^ by setting Vj Â— Â—00 for all j = l,...,n. Let s Â— min {i Â€ {z'l, . . . , ip} : rriij^ = I}. For all j = 1, . . . , n such that there exists an index i E {ii,.. .,ip} satisfying rriij = I, change the value of Vj to a,j + flsji Â— Oiji, where i G {I'l, . . . ,ip}. Note that the value of Vj does not depend on the choice of an element i of the set {ii, . . . ,ip} n {i = 1, . . . ,m : rriij = /}. STEP 4 (Extend horizontal markings vertically and create column vectors u' e RÂ™**^): Define the column vector u' G RÂ™^^ as follows: n = A ^0
PAGE 77
71 Note that u adopts a finite value for alH = 1, . . . , m, since ^4 is a finite matrix and Vj^ is finite. Increase the values of all zero entries of M to / whenever a,j = u[. STEP 5 (Final horizontal extension of markings and creation of final vectors v' 6 RÂ™^^): For all = 1, . . . , n such that is finite, we equate Vj with f j. Otherwise, we define m vj as /\ {aij Â— u). Any remaining zero entry m,j is set to / if aij = Uij + uj. Example. Let us first illustrate how Algorithm 3 works with a simple example. Suppose A G R^^'* is the following matrix: /8 6 3 4" A= 5 3 3 4 \6 6 7 8, A matrix M e R^^"* is initialized as O3X4. Then we execute STEP 1 through STEP 5 for the parameter / = 1. Since 0,3 au = 1 for alH = 1,2, 3, we set (Â«i, . . . , ip) = (1,2, 3), j\ = 3, and 72 = 4. The matrix M G R^^"* now has form /O 0 1 1\ M= 0 0 1 1 . \0 0 1 1 / 3 3 Considering /\ (ozi Â— Cja) and /\ (cji Â— Oja), we recognize that 2=1 z=l A (Â«H a,3) = (8 3) A (5 3) A (6 7) = 1 = asi 033 , 5=1 3 A (Â«^2 "23) = (6 3) A (3 3) A (6 7) = 1 = 032 033 Â• 2=1 Therefore ^0011 M = I 0 0 1 1 .1111
PAGE 78
72 Initially, all entries of the vector are set to oo. Since mu = 1, the parameter 5 = 1, When executing STEP 3, these entries adopt the following values: f J = 031 + Â«13 033 = 6 + 3 7 = 2 . ^2 = 032 + ai3 Â— 033 = 6 + 3 Â— 7 = 2 . ^3 = 0^3 + ai3 aÂ«3 = = 3 Vi = 1,2,3. vl = ai4 + ai3 ai3 = (a,4 0^3) + ai3 = 1 + 3 = 4 Vi = 1,2,3 . STEP 4 computes a vector 6 Rsxi as follows: 4 t}= /\ (ai,t;]) =(82)A(62)A(33)A(44) = 0, 4 1 4 1 4 4 2 = A (o2jt^j) =(53)A(32)A(33)A(44) = 0, ul= f\ (a3j J;) = (6 2) A (6 2) A (7 3) A (8 4) = 4 . Since the vector is finite, we have = in STEP 5. Now the parameter / changes to 2, since there are still some zero entries of M left. These zeros of M are eliminated in STEP 1, since an Â— a\2 = 021 Â— 022 = 2 and mil = "ii2 = "^21 = "^22 = 0. Thus, /2 2 1 1 M = I 2 2 1 1 Vl 1 1 1, Note that STEP 2 through STEP 5 effect no change in the matrix M, since M > O3X4. STEP 3 and STEP 4 create the following vectors e R^""* and u^* G R^""^: v2 = (8 6 Â—00 Â—00 ) ,
PAGE 79
73 STEP 5 produces a row vector e R^^"*, where vf = vf for i = 1,2. The entries vj and vl are computed as follows: vl = /\ (a.3 uj) = (3 0) A (3 (3)) A (7 (2)) = 3 , 1=1 3 vi = A (Â«.4 = (4 0) A (4 (3)) A (8 (2)) = 4 . 1=1 Since all entries of M are strictly greater than 0, Algorithm 3 finishes yielding the following result: A = (u'Mv^)V(u=Mv=) ^0^ 0 I 0(2 2 3 4) V 3 I 0(8 6 3 4) ^8 6 3 4' 5 3 3 4 ,6 6 7 8 Theorem 4.6.5. Suppose A > xBy, where A G R"*^". x G R'^^S and y G R'^". Let u, V, V denote the vectors u', v^, constructed from the matrix A in Algorithm 3 for an arbitrary, but fixed parameter I. Let j\ , j2 be as in Algorithm 3. If I and J denote the sets / = {i = 1, . . . , m : Ui + Vj^ = Xi + yj^ = aij^ and Ui + vj^ = Xi + yj^ = a.j^} f J = < j = 1, . . . , n : aij aij^ = ^ (a^i a^yj for some i G / I 2=1 then O'ij > Ui + Vj = Ui + Vj > Xi + yj V I = 1 , . . . , m; W j E J . In particular, if xi + t/j = aij for some i G {1, . . . , m} and some j G J, then ui + vj = aij.
PAGE 80
74 Proof. Note that / C {ij, . . . , ip} since a^jj aij^ = Vj^ vj^ for alH 6 /. Let us first show that Qij > Ui + vj > x, + yj for all z G I and for all j e J. By Lemma 4.6.3, it suffices to show that for all j Â€ J there exists an index ui^ + vj Â— ai^j > xi^ + yj. By the definition of J, each index j 6 J entails an index ij G / such that mi^j is set to / in STEP 2 of Algorithm 3. If s denotes min{i G {ii, Â• ,Â«p} : rnij^ Â— /}, then Vj = ai^j + Qsj^ a,jj^ for all j G J. For an arbitrary, but fixed element j E J the inequalities vy < ai^y + agj^ Â— ai^j^ hold for all / = 1, . . . ,n. Therefore, we are able to deduce the following equalities: n = (a, .J Vj) + = Oiji VjGJ. Since A > x El y, the conditions of Lemma 4.6.3 are satisfied, yielding that u, + vj > xi + yj for all? G / and for all j G J. Provided that w,' + Vj^ = aij^ or u,' + = aij^ for all i' G {1, ...,Â«} \ / can be shown, an application of Lemma 4.6.4 concludes the proof of the theorem, since Vj = vj for all j Â£ J and yl > u El v by the definition of u. Let us fix an arbitrary row index n i' ^ I. By the definition of u,as /\ (a,/j Â— vj), there exists a column index j' such j=i that Ui' Â— ui'j' Â— Vj'. Since u,/ is finite, Vj' is finite and equals vj'. Therefore, j' either equals j2 or is an element of the set J' which we define as the set of all column indices m j such that Ufj a^jj = /\ {a^j Â— a^j^) for some z G {^i, . . . , ip}. Note that J' includes 2 = 1 jl. The case j' = j2 leaves nothing to show. Assuming that / G J', we finish the proof of the theorem as follows. We use the facts that A > u El v and that a,jj = + vj^ f
PAGE 81
75 for all i e {h, . . .,ip}. m ai'j' Gi'j^ > f\ {a^ji a^j^) = a,j' a,jj for some z Â€ {ii, . . . , ^p} 2=1 ^ Ui> + vji ai'j^ > {ui + Vji) {ui + Vj^) for some i E {h,.. .,ip} Ui' + Uji > Cj/ji Ui> + {jjj = ai'j^ . Theorem 4.6.6. Algorithm 3 is a polynomial time algorithm. Proof. Let us analyze STEP 1 through STEP 5 for A G R'"''" and for fixed / Â€ N. Since there are "^"^"^^ different tuples of column indices and j2, m lii^iH subtractions are needed in order to compute the differences aij^ Â— aij^ for alH = 1, . . . , m and for all ji,j2 = I,,". For fixed ji and j2, the comparison of the m real numbers a,jj Uij^ takes "^^^"^^ operations. Hence, these comparison operations are of computational complexity 0{\ Â• ri^m^). All other operations of STEP 1 and the operations of STEP 2 through STEP 5 have smaller computational complexity. Since the number of zero entries of M Â€ R"*^" decreases when performing STEP 1 through STEP 5, Algorithm 3 is a polynomial time algorithm. Theorem 4.6.7. Algorithm 3 produces the desired results. Proof. Suppose that Algorithm 3 is applied to a matrix A 6 R"*^, yielding as an output an integer k, column vectors u' e RÂ™^^ and row vectors v' Â€ R^^** for all / = 1, . . . , k, as well as a matrix M E R"*^". Since the final matrix M belongs to the set {1, . . . , /j}"*^", it suffices to show that the following properties hold for all i, j, and /: 1. u\ \v'j < Uij.
PAGE 82
76 2. rriij = I ^ u\ { Vj > aij. n The first property immediately follows from the fact that u equals /\ (a,^ ) , since n A ^z) < (^a z=l Suppose that m,j = / after completion of STEP 3. This assumption implies that i e {ii,...,ip} and Vj = = a,j + agj^ a,ji. In this case, the following inequalities show that property 2 holds: n "! + "i = A 2=1 n > /\ (Â«i2 (a,2 + dsji aiji)) + '^j 2 = 1 = ) + (a.j + Â«sii Â«tii ) Clearly u'+ v'j = aij, if the entry m,j of the matrix M Â€ R"*^" has been set to / in STEP 4 or in STEP 5. Theorem 4.6.8. Given a matrix A G R"*^" of rank 1, Algorithm 3 determines a column vector G fjmxi ^ vector E R^^", 5mc/i that Proof. By Li's Theorem, each row vector a{i) of A, where i = 2, ...,n, is linearly dependent on the first row vector a(l). Hence, there exist such that a(i) = + a(l), which implies that the following equalities hold: an o,\2 = A, + ail (Ai + 0,2) = a,i a,2 V z = 2, . . . , n .
PAGE 83
77 Hence, STEP 1 of Algorithm 3 results in a matrix M 6 R"*^" of the following form: M = /I 1 0 0 110 0 0\ 0 yi 1 0 0 0/ Since the differences a^j a^i , where z ranges from 1 to m, are all equal for arbitrary, but fixed j e {I,. . . ,n}, all entries of M are set to 1 in STEP 2. STEP 3 yields a row vector = a(l). In STEP 4, we compute a column vector as follows: = an Oil Â• Hence, = a[l] an. STEP 5 equates with and the algorithm finishes yielding A = u^Elv'. Example. The following example shows that although Algorithm 3 always determines a weak decomposition of a given matrix A into outer products of column vectors and row vectors, the number of these outer products is not always minimal. In other words, there are instances when the parameter r determined by Algorithm 3 exceeds the rank of the matrix A. Let A 6 R^^^ be the following matrix. n n
PAGE 84
78 /2 0 0 1 4\ 5 3 10 8 3 5 4 2 7 \3 1026/ V /o \ 1 3 \v /4\ 0 2 V2/ 0 4 Vo / El(5 0 0 1 4) V 0(5 3 1 0 8) E](3 1 0 2 0) Since A is weakly decomposable into three vector pairs, the rank of rank(v4) < 3. Using Algorithm 2, we can show that rank(yl) = 3. However, Algorithm 3 yields k = i and computes the following column vectors u\ . . . , and row vectors v\ . . . , such 4 that A = V u' El v'. /o\ /7\ /o\ 3 0 ,u3 = 1 5 2 4 v) \0j ,U4 = /O \ 1 3 v^ = (2 0 2 3 2), = (5 3 1 0 8 ) , v3 = (2 0 0 2 3), v4 = (2 0 0 1 4). Computational Results. We have coded Algorithm 3 in MATLAB. We randomly generated 100 matrices of size 16 x 16 having integer values ranging from 0 to 255 and rank
PAGE 85
79 r for r = 3,4,5. We succeeded in finding a weak decomposition of the given matrices into r vector pairs for the following number of matrices: 1. rank 3: 91%. 2. rank 4: 88%. 3. rank 5: 74%. The average CPU time and the average amout of flops required for the rank decomposition depend on the rank of the given matrix A: 1. rank 3: 13.45 sees, 75885 flops. 2. rank 4: 18.02 sees, 110950 flops. 3. rank 5: 24.03 sees, 146710 flops. Note that the intended use of Algorithm 3 is the problem of weakly decomposing a morphological m x n template into r outer products of column templates and row templates. Savings in computation time can only be achieved if r < ^i^.
PAGE 86
CHAPTER 5 THE SVD AND THE MED 5.1 Introduction In the previous section, we discussed rankbased decompositions of matrices in the minimax algebra domain, intended for the use of morphological template decomposition. Potentially, rankbased matrix decompositions both in the linear and in the nonlinear domain could also lead to image compression techniques provided that the rank of the matrix corresponding to the image is relatively small compared to the image size. Reallife images have prohibitively large rank in linear algebra as well as in minimax algebra. However, image compression using outer product expansions does not require an exact decomposition of the corresponding matrix but merely an approximation of this matrix by a matrix of relatively small rank in terms of a suitable metric. In the linear domain, the singular value decomposition (SVD) provides the best approximation of a matrix A by a rank k matrix in the least squares sense for each k < rank(>i) [28]. Barring illconditioning of A, excellent approximations of A can be achieved even if k is significantly smaller than rank(A), making the SVD an ideal candidate for image compression purposes. Furthermore, singular value decomposition methods have proven useful in image analysis, image representation, image storage, and image restoration [1, 2]. The definition of the rank of matrices in minimax algebra in terms of outer product expansions enables us to determine a minimax algebra analogue to the SVD, which we call minimax eigenvector decomposition, or briefly MED [51]. After presenting this new 80
PAGE 87
81 image transform, we conclude with a brief comparison of properties of the MED and the SVD. 5.2 The Singular Value Decomposition The singular value decomposition (SVD) of a realvalued m x n (m > n) matrix A allows for the representation of A as a matrix product A = USV'^, (51) where U and V are orthogonal matrices of dimension m x m and n x n, respectively, denotes the transpose of V, and 0 0 Â• Â• 0 0 02 0 Â• Â• 0 0 0 0 0 Â• Â• 0 0 0 0 Â• Â• 0 0 \o 0 0 Â• Â• 0 0 / is a diagonal matrix whose diagonal entries are ordered in the form (7 > Â• Â• Â• > CTnSpecifically, the cr, 's are the square roots of the eigenvalues of the matrix A'^A and are called the singular values of A. The columns of U = (u\ . . . , u*") are the eigenvectors of AA^ and the columns of V = (v\ . . . , v") are the eigenvectors of A^A. The vectors u\ . . . , u"* are known as the left singular vectors and the vectors v\ ... , v" are known as the right singular vectors. A well known application of the SVD in image processing is in image compression. Observe that Eq. (51) can be rewritten as r A=;^^.u'(v')'^, (53) 1=1
PAGE 88
82 where r denotes the rank of S and (v') denotes the transpose of the vector Vj. Image compression can be achieved due to the fact that (5^) where /; < r. Thus in order to store the whole m x n image, all we need to store is the set {cr,, u', v'l^^j. Suppose that 1 1^"! I2 symbolizes the familiar 2norm of an arbitrary realvalued matrix X. It is wellknown that ^2 equals the largest singular value of the matrix X. The following theorem shows that a SVD of a matrix A provides the best approximation of A with respect to the metric induced by this norm. Theorem 5.2.1. Let A^ UEV* be a SVD of a matrix A Â€ R'"''". Setting k Ak = Y,
PAGE 89
83 methods for finding eigenvalues are all iterative. The amount of work required in each iteration is often prohibitively high unless initially A (or AA*) is of special form. One of the best methods available for computing the singular values of a matrix is as follows. Observe first that if A has singular value decomposition A = U EV^ and B is a matrix of form B = HAP^, where H is an orthogonal m x m matrix and P is an orthogonal nxn matrix, then B has singular value decomposition {HU)E{PV) . Thus the problem of finding the singular values of A can be simplified by applying orthogonal transformations to A in order to obtain a simpler matrix B with the same singular values. Golub and Kahan have shown that A can be reduced to upper bidiagonal form using Householder transformations and an excellent exposition of this method is given in [28]. This brief sketch of the methodologies for finding singular values should provide an appreciation of the complexity and intricacies of the eigenvalue problem. An additional problem concerns the accuracy of the eigenvalues one obtains by using any of the current methods. Since these iterative matrix solvers require many multiplications and divisions for large matrices, the error cr, Â— a\\ between the computed and actual singular values may be quite large. In contrast, the computation of eigenvectors for minimax eigenvector decomposition does not involve these types of complexities and errors. 5.3 The Minimax Eigenvector Decomposition (MED) An element A G RÂ±oo is a minimax eigenvalue of A with corresponding eigenvector X Â€ (RÂ±oo)", if AHx = A + X. Here A + x=(A + xi, ...,A + a;Â„), where x = (xi, . . . ,a;Â„).
PAGE 90
84 Theorem 5.3.1. [14] If A e RÂ±(^", then there exists a matrix U G R'^'^ such that each column of A is an eigenvector of U corresponding to eigenvalue 0, namely U = A^A*. In other words UUiA = A. Proof. It suffices to show each of the inequalities UEBA > A and U M A < A. An bound for {U M A)^j can be obtained as follows: {U M A)^j = y (uih + ahj) h=i h=i n V h=l t=l k=l m > /\ {aik +' {aik)*) + aij k=l > aij = l,...,m; V; = l,...,n . The last inequality holds because for any a G RÂ±oo the sum a +' (a)* equals 0 for finite a and oo else. Similarly, a\a* = a*\a={ " . This fact and the first Â— oo else associative inequality in Theorem 3.2.2 lead to a proof oi U MA < A similar to the
PAGE 91
85 above proof of the inequality UMA> A. <{Am{A*MA))^^ = A h. n = A < h=l L fc=l n dih +' V {{akh)* + akj) n h=l L k^l n jfc=l < ttij Vi = l,...,m; Vj = l,...,n . Therefore, C/ 13 A = A and we have certain parallels to the SVD of linear algebra. Specifically, the MED of A has form A = UmSMV\ (56) where V = A*, oo S = oo Â•72 Â— OO Â— oo Â— oo Â— oo oo oo oo oo oo Â— oo Â— oo Â— oo Â— oo oo \ oo oo (Tn oo Â— oo oo Â— oo , and (Ti = (T2 = Â• Â• Â• =
PAGE 92
86 to form the conjugate matrix A* and compute the product U = AEiA*. In analogy to Eq. (53), it follows from Eq. (56) that n 1=1 (57) Since (v*)' denotes the ith row vector a' of A and cr, = 0 for Â« = 1,2, . . . ,n, we actually have the simpler formulation n A=\/ {u'M a') . (58) i=l Since the u"s are the columns of the matrix U = A\S A*, the determination of the components of the u"s requires only additions and the logic operation A of minimization. The components of the row vectors a' of A are given and the computation of the eigenimage u' El a' requires only additions. Thus the determination of eigenvectors and eigenimages required in the minimax decomposition of a matrix is for all practical purposes trivial. The following example serves to demonstrate this observation. Suppose 9 5 Then A* ^ \ 7 4 1 \8 7 0 and, therefore. ^0 15 8^ U = A^A* = \ 6 0 4 ,2 7 0 5.4 Comparison of the MED and the SVD In comparison to the SVD, minimax decomposition of a matrix is trivial. However, what is often desired is not the complete set of outer product expansions but only a small
PAGE 93
87 subset of these with the property that this subset provides for the best approximation Ak of the original matrix A. In case of the SVD, once the singular values and singular vectors for the series expansion of Eq. (53) have been determined (the difficult part), it is trivial to obtain the best approximating subset of k outer product expansions for a given k
PAGE 94
88 2. Suppose the first j (i > 1) eigenvector pairs have been found. Find two vectors and a* such that (u*Ela") V V(u Ela^) 1=1 J 1=1 = y d, i=j+l (u' 0 a') V V (u' 0 a') /=i ,\/(u 0a')) . J 1=1 Rename this pair as u'"'"^ and a'"''^. 3. Set i = i + 1 and continue with step (2.) and step (3.) until j = k. k An approximation of A with respect to the norm di is given by Ak = V 0a'). Note that if A 6 Z"*^", then all necessary computations can be performed within Z. Hence, this algorithm guarantees numerical stability for the case A Â€ Z"*^". Example. If we apply Algorithm 4 to the 5 x 5 matrix /39 7 0 94 84\ 62 0 58 11 3 A = 73 19 30 2 0 9 75 8 20 29 \58 27 18 82 4 / then using the first three eigenvector pairs we obtain the reconstructed matrix /39 7 0 94 84 \ ^3 = V (u' 0 a') = Â»=i 62 0 58 11 3 34 0 30 2 0 9 0 5 20 10 \58 27 18 82 4 / Using the complete set of five eigenvector pairs we obtain, of course, A^ Â— A. To contrast this with the SVD, the series expansion A3 = ^ cr,u'(v')^ yields 1=1 /36.813789 8.839371 2.534787 103.915436 70.123454 \ 69.110519 1.802377 43.164482 15.181345 7.691323 A3 = 63.394817 19.603844 41.031982 10.571262 7.335522 10.670452 75.780029 3.457618 22.650938 23.932810 V60.955730 23.081055 25.672878 61.892658 33.012875 /
PAGE 95
89 while As /39.000019 6.999990 0.000003 94.000031 84.000000 \ 62.000015 0.000008 58.000008 11.000007 2.999976 73.000015 18.999996 30.000013 1.999997 0.000021 9.000000 75.000031 8.000004 20.000019 29.000017 \58.000008 26.999996 18.000004 82.000008 4.000021 / Note that although in theory A5 = A for the SVD as well, in actual computation is only approximately equal to A due to numerical errors introduced by finite digital computers when using multiplications and divisions. Example. As an additional example, we compare the MED and SVD transform on an actual image. Figure 9 shows the 12 x 12 input image, while Figures 10 and 11 show the eigenvector expansions for three and seven eigenvector pairs of the SVD and the MED, respectively. Figure 9 The 12 x 12 source image A. Figure 10 Image reconstruction using the SVD transform; the image on the left corresponds to A3 and the image on the right to Aj.
PAGE 96
90 Figure 1 1 . Image reconstruction using the MED transform; the image on the left corresponds to and the image on the right to Aj. Examples like these indicate that image reconstruction using the MED works well for "boxy" images of relatively small size. However, if A G R"Â»xÂ« corresponds to an m X n real life image, the distances di{A,Ak) are prohibitively large for all A; < n. Moreover, Corollary 4.6.2 implies that the matrix Ak does not necessarily represent the best approximation with respect to the metric di of A among all matrices of rank k. We dealt with these problems concerning the image reconstruction of reallife images using the MED by employing the technique of divide and conquer. Let A 6 R"*^'*. If p divides m and q divides n, then A has the following block structure, where A* eRp A = / A' A= ... Ai\ ^^(pi)g+i API J For given k
PAGE 97
91 Example. Let us apply this technique to the following realvalued image of size 32 x 32, which corresponds to a matrix A G r32x32 j^iq matrix a[^''^^ can be stored on a digital computer using only a quarter of the entries of the matrix A Â€ r32x32^ since the representation of the matrix A\ = u* M (v*) G R^^^ requires only 2 Â• 8 = 16 entries compared to the 64 entries of an 8 x 8 matrix. Similarly, forming Ai^''*'^ yields an image compression rate of 50 %, since the matrices A\ G R^^^ have representation in terms of 2 outer products of vector pairs of length 8. Figure 12. The 128 x 128 source image A. Figure 13. Image reconstruction using the SVD transform; the image on the left corresponds to Aie and the image on the right to A32.
PAGE 98
92 Figure 14. Image reconstruction using the MED transform; the image on the left corresponds to about one fourth and the image on the right to half of the total information.
PAGE 99
CHAPTER 6 MORPHOLOGICAL DECOMPOSITION AND APPROXIMATION 6.1 Introduction The previous chapter describes the MED transform and it's use for determining the best approximation of a matrix with respect to the metric di. By Corollary 4.6.2, the outer product expansions given by the MED transform do not necessarily yield the best approximations of the original matrix. In this chapter, we are presenting an integer programming approach to matrix approximation problems using di as a distance. By the definition of a metric, matrix approximation in terms of outer product expansion results in an exact representation of the original matrix, if the distance between the norm and the approximating outer product expansion is 0. Due to the natural correspondences between matrices and images as well as between matrices and templates, solving matrix approximation problems could have dual use for image and template applications. However, since the size of images would lead to inmiense optimization tasks, this approach is only intended for the purpose of template decomposition and approximation. Using this approach, we are able to tackle the problem of weakly decomposing arbitrary given templates with respect to the operation El in its most general form. We provide exact definitions below. 6.2 Morphological Decomposition and Approximation Problems The decomposition problems, we are primarily interested in, are matrix and template decompositions with respect to the operation El . Furthermore, we assume that the 93
PAGE 100
94 templates which are to be factorized are elements of (R?oo) where X is a finite subset of 1? (without loss of generality X contains z := (0, 0). Linear mixed integer programming is suited to solve a very general form of these problems. The General Form of the Decomposition Problem. Let A e R'"^". Assume that A is decomposable as follows for some natural number h: A = BUC, C e R*^''" . Suppose a set Y is defined as follows: Y = {yi, Â• Â• ,ym*fc+fc*n} = {^ll,,'>lfc,,^ml,^mik>Cii,...,CiÂ„,...,C;tl,,CjfcÂ„} Then there are linear functions fjj {i = 1, . . . ,m; j Â— 1, . . . ,n; / = 1, . . . , A;) such that for each value of i and j the following inequalities hold: fijiVl, ,ymk+kn) < aij fijiVliÂ• Â•,ymk+kn) < (ij) fij{y\^ ^ymk+kn) < Â«Â«> Â• Moreover, for each i, j there exists an index /'^ e {l,...,k} such that: fij (yi, Â• Â• Â• , ymk+kn) oij Obviously the linear functions fj are simply derived from the additive maximum B 0 C: fljiyii5 ymk+kn) = bii + Cij .
PAGE 101
95 A similar result holds for any weak El decomposition of an invariant template t Â€ (R?"oo) with finite support into invariant templates t%...,t*^S...,t*^"^''"\...,t^" G (Roo) Â• The support of each template must be given in advance. Suppose the following definitions are given as well: X = Xj {xi,...,Xm} = S{t) Y = {yu . . .,ym'} = t'{S{t')) U . . . U t^" ("^ (*'")) Â• For each i = 1, . . . ,m, there exists a natural number h and linear functions // for all I = 1, . . . , /, satisfying inequalities of the following form: //(2/i,,J/m') < tx(x,) /f(yi,,ym') < tx(Xi) [i] fHyxi'iVm') < tx(x.) . Furthermore, there must be at least one ki Â€ {1, . . . , Z^} which allows us to write one of these inequalities as an equality: ftiyi^^ym') = tx(x,) . The problem of finding such a decomposition could be solved as suggested by Takriti and Gader who are particularly interested in weak decompositions of invariant templates into 3x3 rectangular templates [24]. However, we have taken a more general approach. Moreover, we are also able to find the best approximation of a template with respect to a certain norm.
PAGE 102
96 We introduced the reader to the concepts of matrix and template decompositions with respect to El (we also call them morphological decompositions). Now it is time to give exact definitions of the problems we are concerned with: i. Morphological matrix decomposition problem: Given a matrix A E R"*^" and k Â£ N, determine if A is weakly decomposable into k vector pairs. If possible, find B G R'"^'= and C e R^^" such that BMC = A. ii. Morphological template decomposition problem: Let t G (Roo) ^ ^ invariant template, (fci, ^2, . . . , /^n) an strictly increasing sequence of natural numbers and (X^X'^, . . . ,X^") a sequence of subsets of X. Find out if there exist invariant templates t\ t^ . . . , t*" G (R^oo) such that: a. For alH = 1, ...,/?Â„ the set X' contains the support of t^, where z = (0, 0). b. The template t factorizes as follows into a product of the templates , , . . . , t*^" : t = ^t' El . . . E t*^^) V . . . V ^t*="'+' El ... El t*") . Determine the templates t\ t^, . . . , t*^" in this case. These decomposition problems are closely related to the following approximation problems with respect to the distance di. I. Morphological matrix approximation problem: Given a matrix A G R"***" and A; G N, find B G R'"^'^ and C G R^""" such that di{BMC,A) < di{B'MC',A) V5' G R'"''^ C G R^""" .
PAGE 103
97 II. Morphological template approximation problem: Suppose t Â€ (R?oo) ^ invariant template, {k\,k2,. . . ,kn) an strictly increasing sequence of natural numbers and (X^X'^, . . . ,X^'") a sequence of subsets of X. Define T to be the set of all sequences s^, . . . , s*") of invariant templates which have the following properties: 5(s' ) c X' Vi = 1,...,A;Â„ . Find a sequence of templates (t^,t^, . . . ,t*^") e T which satisfies the inequality (fi ( El . . . El 1*=^ V...V t'="^+^E ... E!t*"],t) < ^ El . . . El s*^] V . . . V M ... El s*"] , t) v(sS...,s*^") ET Since we have the equivalences di{A, D) = 0 A Â— D and c?i(t, r) = 0 <^ t = r for all matrices D Â€ R"*'^" and all invariant templates r Â€ (R^oo) Â» *he solutions of the approximation problems lead to solutions of the decomposition problems as well. The existence part of the morphological decomposition problems defined above can be answered depending on the distance of the best approximation to the original matrix, and the original template, respectively. If this distance in terms of the metric di equals 0, then the approximation problems and the decomposition problems have the same solutions. 6.3 Integer Programming Formulations Consider the following slight variation of the matrix approximation problem I described in the last section.
PAGE 104
98 r. Suppose that A e R"*^" and G N are given. Define V as the set of all matrices D' = B'm C such that B' g R"'''^ C G R^""" and D' < A. Find matrices B G R'"''* and C G R^""" which satisfy BMC < A md di{A,BMC) < di{A, D') for all V D' G P'. Since we are restricting ourselves to the set we will not find the best approximation overall of the matrix A G R"*^" in terms of an additive maximum BMC. However, we are still able to find out if A can be decomposed exactly as B M C for some B G R"*^* and C G R*^^". In this case, problem I' yields the same solution as the problems I and i. m n For all D' eV the distance di{A, D') amounts to nothing else than ^ ^ a,j 1=1 3=1 d^j . because D' < A. Therefore, it suffices to find matrices B G R"*^* and C G R^^" m n such that D = BMC < A and X) 'S maximal. Â»=i j=i Now let us return to the system of inequalities (ij) given in the first section of this chapter. Notice that for each group of inequalities depending on i, j only one equality must hold. Takriti and Gader call this fact "the group equality restriction" [24]. We introduce variables r'^ indicating which equality constraints in (ij) are fulfilled and deduce the following formulation: k m n ""'J * (^'^ + ""'j) 1=1 i^l J=l subject to rlj = 0 or 1 k bil + cij < aij V Â« = l,...,m, j = I = l,...,k .
PAGE 105
99 Equivalently, this maximization problem (I') can be formulated as the linear mixed 01 integer programming problem given below. k m n /=i i=i >=i subject to a * <
PAGE 106
100 all / = Let us return to problem I. The best approximation with respect to the metric di of a given matrix A G R"*^" in terms of a matrix product BMC, where B Â€ R'"'^*, C e f^kxn^ is requested. Matrices B* and C*, which solve this problem, can be found by solving the following minimization problem. The optimal objective value will equal it Qij V {bii + cij) di{A,B*EliC*) =EE i=l j=l 1=1 ,=1 j=l \ 1=1 subject to p\j > aij kl cij qij > hi + c/j a,j qij > 0 rjj = 0 or 1 E r'= 1 V i = 1, . . . ,n, J = 1, . . . ,m, / = We would like to point out some helpful observations to the reader: Notice that each variable p\j is greater than or equal to 6,/ c/j) , the positive part of aij {bn + cij). The variables r\select an index /' G {1,...,A;} such that (a,j 6,/' c/zj)"^ is minimal, or equivalently bui + ciij is maximal whenever 6,7+c/j < a,j. Consequently,
PAGE 107
101 for each i,j the optimal solution satisfies k / k \ + * Pij = ( Â«o V (^Â«' + ^o) ) Â• Â• The variables qij exceed the values 0 and 6,7 + c/j Â— a,j for all / = 1, . . . , A;. Hence, at optimality A linear 01 MIP formulation can be given as follows: m n / k :=1 j=l \ /=l subject to Pij > aij hii cij fi * s\j Qij > hi + c/j aij qij > 0 5I, = 0 or 1 V i = 1,. . . ,m, j 1,. . . ,n, I = l,...,k . The variables sj^ take on the role of the variables r'^ from the quadratic mixed 01 approach. The number fx should be at least as large as V a,j Â— (bn + c/j). Similar 01 MIP formulations exist for the corresponding template approximation problem II. These formulations are based on systems of inequalities of the form [i].
PAGE 108
102 Mixed 01 Integer programming can also be applied when the maximum metric d2 instead of di is used. Example. Let A; = 2 and A Â€ R^*^^ be the following matrix: /7 12 11 A = (9 18 9 \5 14 7 By the Example given in Section 3.3, the matrix A can be written as A = BMC, where B e R^^" and C Â€ R^^^ Therefore, the 01 MIP program below will yield the optimal 3 3 objective value ^* = Z) Â— 92: Â«=i j=i 2 3 3 l^l i=i j=i subject to 18 * rlj < p\j < 18 * r!_, r[j = 0 or 1 p\j < bil + cij h + cij < aij V i = 1,2,3; j = 1,2,3; / = 1,2 . In most cases, we only have to deal with images which only assume nonnegative values, k k i.e. nonnegative matrices. Therefore, we can replace ^ r' = 1 by r< 1, which /=i /=i will replace negative entries in B El C by zeros and will thus improve the approximation of A.
PAGE 109
103 The sizes of our mixed 01 linear programming problems increase with the size of the matrices and templates which are to be decomposed. For example, the morphological matrix approximation F for an m x n matrix and k vector pairs has k{2mn + m + n) variables and {ik + l)nm constraints (cf. the preceding example). We tested our linear 01 MIP formulations on some matrices of sizes 8x8 and 16 X 16 by using the commercially available optimization package LINDO [53]. In spite of very limited CPUtimes, we were always able to obtain good approximations of the given matrix (i.e. an objective function value within 10% of the exact solution). A detailed description of two specific examples is included at the end of this chapter. In the course of the last few decades, many efficient 01 mixed linear algorithms have been developed and implemented to solve larger problems [12, 26, 45]. More recently, interior point approaches have been successfully used to solve very large 01 mixed integer problems optimally or approximately. The success of these approaches is due to very efficient implementations of interior point methods for solving very large scale linear programs [34]. We concentrate our current efforts on solving computer vision related mixed 01 integer problems with interior point methods. Extensive computational results with interior point methods will be presented in a forthcoming article. Some Preliminary Results: Decomposition of a Matrix into 3 Vector Pairs with LINDO. Let the following matrix A G R^^^ be given. This matrix is known to be
PAGE 110
104 weakly El decomposable into 3 vector pairs. /19 15 17 9 20 19 21 24 13 23 18 14 18 7 18 21 17 21 13 24 22 18 22 20 31 25 21 21 9 15 24 20 20 8 17 \28 24 24 16 27 15 14 15\ 18 21 22 13 15 16 19 18 19 26 22 21 16 17 10 15 16 12 22 20 17/ The corresponding 01 MIP formulation for problem I' is characterized by the following data (cf. Example 9): Â• number of 01 variables: 3 * 64 = 192 Â• total number of variables: 3*(2*64 + 2*8) = 432 Â• number of constraints: (3 * 3 + 1) * 64 = 640 Â• exact solution: 1182 Â• approximate integer solution (5 min CPU time on VAX8600): 1073 Â• number of pivots: 2451 Finally we present some data and results of an application of the MIP algorithm to a matrix B G Ri^xie matrix B also factors into a El product of 3 vector pairs. Â• number of 01 variables: 3 * 16^ = 768 Â• total number of variables: 3 * (2 * 16^ + 2 * 16) = 1632 Â• number of constraints: (3*3 + 1)* 16*^ = 2560 Â• exact solution: 4256 Â• relaxed linear programming solution ( 15 min CPU time on VAX8600): 4630 Â• number of pivots: 2116
PAGE 111
CHAPTER 7 COMPLEXITY OF MORPHOLOGICAL DECOMPOSITION PROBLEMS 7.1 Introduction In this chapter, we are going to prove the NPcompleteness results which we addressed in the previous chapters. In particular, we show that the morphological matrix decomposition problem is NPcomplete. 7.2 An NPComplete Problem Theorem 7.2.1. The following problem IT is NPcomplete: INSTANCE: Let A G R'"''" and A; be a positive integer. QUESTION: Is there a collection of vectors u' 6 R'"^* and v' 6 R^^", / = 1, . . . , fc. Proof. (1) n e NP: For any guess of vectors u' 6 R'"^* and v' Â€ R^**", the checking stage can be performed in polynomial time: 1. Forming V (u' 121 v'^ costs k{mn) + {k l){mn) = {2k l)mn operations. such that k 2. Another mn operations are needed for comparing the result with A. 105
PAGE 112
106 (2) Choice of an NPcomplete Problem 11': INSTANCE: Collection C of nonempty subsets of a finite set S and a positive integer QUESTION: Is there a collection B of subsets of S with 5 = ^ such that, for each The following argumentation shows that 11' is indeed NPcomplete: Since IT' represents a subproblem of SET BASIS which is known to be NPcomplete, it is contained in NP. The problem of SET BASIS only differs from 11' by allowing the empty set as an element of C [25]. The problem of SET BASIS can be partitioned into the problem 11' and a similar problem 11". Instances of 11" consist of those collections C of subsets of S which include 0, and a positive integer k < \C\. The problems IT' and n" are polynomially equivalent. Thus, if either one would be contained in P, the problem of SETBASIS would belong to P. (3) Construction of a Transformation / from 11' to 11: Let an instance (C, k) of II' be given, where C is a collection of subsets of a finite set S and A; is a positive integer with k < \C\. The cardinality of C is bounded from above by 2l^l, the cardinality of the power set of S. Hence, C as well as S are finite sets, and they are of the form C = {Ci, . . . ,Cm} and S = {^i, . . . , 5Â„} for some m, n G N. We now construct an instance {A, k) of 11 which is satisfiable if and only if the instance (C, k) of n is satisfiable. The parameter k remains unchanged, and we define the matrix A G R!!*<^" by setting k < \C\. C G there is a subcollection of B whose union is exactly C.
PAGE 113
107 "<=": Let (C, k) be satisfiable, where C is a collection of subsets of S = {si, . . . , 5Â„} and kisa positive integer with k < \C\. Hence, there exists a collection 5 = {Bi,...,Bjt} of subsets of S = {^i, . . . , Sn} such that, for each C, E C, where i = 1, . . . , m, there is a subcollection 5, C B whose union equals C,. We are now going to show that the instance (A, k) of 11 is satisfiable by constructing row vectors Â€ R*^" and column vectors u' Â€ R"***^ for all / = l,...,k, such that it A = y (u'eIv'). 1=1 We define the row vectors v' Â€ R^^", / = l,...,k, as follows: ^^ = lo else V; = l,...,n Introducing the notation L, = {/ 6 {1, . . . , ^} : B/ Â€ Bi} for all i = 1, . . . , m, we define the column vectors u' G R"*^\ / = I,. . .,k, by setting / JO if I eU vj 1 u\ = \ Â•> ' Vz = l,...,m. Â— 1 else k Note that the matrices A and V (u^Elv^) are elements of {0, i}'"'*" c R"Â»xn. /=i it Hence, in order to show that A = V (u' El v') , it suffices to show that /=i {{ij) G M X N : a., = 1} = 6 M x N : \/ (u^Elv^).^. = l , where M denotes from now on the set {1, . . . , m} and N denotes from now on the set {1, . . . , n}. Furthermore, we set K = {1, Â• Â• Â• , k}. The above equation follows from the sequence of equivalences below, where i 6 M
PAGE 114
108 and i Â€ N are arbitrary indices. V(u'0v') =1 k /=i 3 / G K : v'j = l and u = 0 <^ 3 / G L, : Sj e B/ <^ sj G B/ /or some / Â€ <^ Sj Â€ C,4^ Â— 1 . We now conclude the proof of the desired equality A Â— \/ (u^ M v') by showing that ib the entries of the matrix V El v') Â— like the entries of the matrix A Â— only adopt /=i k the values 0 and 1. In other words, it suffices to show that V El v') . e {0, 1} for /=i '^ all i Â€ M and for all j Â€ N. Recall that the sets C, are not empty for alH Â€ M. k V i G M 3 i, G N : 1 = aij, = \/ (u' M v') = 1 V i G M 3 i. G N, /. G K : + = 1 =J> V i G M 3 i, G N, /, G K : u' = 0 =4> V i G M 3 G K : uj+ uj' > 0 V i G N k =^ V i G M : V f 121 v') > 0 V i G N k (n'My^y >0 V i G M, V i G N .
PAGE 115
109 "=j'": Let {A, k) be a satisfiable instance of 11, which means that there exist vectors k g Rmxi gjjj ^/ g Rixn^ / = 1, . . . , A;, such that A = V M v'). We need to show /=i that the corresponding instance instance (C, k) of 11' is satisfiable. We will construct a collection B of k subsets of S such that, for each Ci 6 C, i = 1, . . . , m, there is a subcollection Bi C B whose union is exactly Ci. For all / e K, we define B/ as the set of all s^' Â€ S such that vj, is maximal among all Vj, where G N and denote the collection of these sets by B, i.e. 5 = {BÂ„...,Bfc}. If L, C K, z = 1, . . . ,m, denotes the set of all / G K, such that u'+ vj = 0 for some j G N, then we define By the definition of A G R , Ci = {sj G S : aij = 1}. Therefore, /et;
PAGE 116
110 (4) Proof that / is a Polynomial Transformation: By the construction in (3), / maps an instance (C, k) of 11', where C = {Ci, . . . , Cm} is a collection of m subsets of S = {si, . . . , 5Â„}, to an instance (A, k) of 11. Since the parameter k does not change, the only operations involved in the computation of / are the ones arising from the computation of A G R"*^". For each entry a,j of A, where i G M and j e N, the element Sj of S has to be tested if it belongs to the set C^. This test involves less than n operations. Thus, the computation of / can be performed in less than mn^ operations which implies that / is a polynomial transformation. 7.3 Extensions Finally, we provide a few easily derivable, but important corollaries to the above theorem. Corollary 7.3.1. The problem of finding the rank of a rectangular morphological template is NPcomplete. Corollary 7.3.2. The problem of finding a rankbased decomposition for rectangular morphological templates is NPcomplete. Corollary 7.3.3. The general morphological template decomposition problem is NPcomplete. Concluding Remarks. We attempted to provide a solid theoretical background for the complexity of all morphological template decomposition problems. We have shown that a particular subproblem 11 of the general morphological template decomposition problem is NPcomplete. Hence, global optimization approaches to solve the morphological template decomposition problem can not employ linear programming techniques. Attempts to
PAGE 117
Ill exploit special structures of certain morphological template decomposition problems will continue to play a predominant role in future efforts to efficiently solve morphological template decomposition problems.
PAGE 118
CHAPTER 8 CONCLUSION AND SUGGESTIONS FOR FURTHER RESEARCH This dissertation represents the first account of matrix decomposition methods in terms of outer product expansions in minimax algebra beyond the decomposition of separable matrices. Due to the correspondences between matrices and templates and between matrices and images, these methodologies lead to applications in image processing. We provided a solid mathematical basis for matrix decomposition by establishing a new theory of rank within the framework of minimax algebra. The resulting concept of rank resembles the concept of rank known from linear algebra in the following respect: A matrix A has rank n if and only if n is the minimal number of outer products which allow for a representation of A. Given this background, we established the following results: 1) We related our definition of matrix rank in minimax algebra to the one given by CuninghameGreen. We showed in particular that the rank of a realvalued matrix A is bounded from below by the rank of A according to CuninghameGreen and bounded from above by the maximal number t of linearly independent row vectors of A as well as by the maximal number of linearly independent column vectors of A. 2) We generalized Li's Theorem which provides for the strong decomposition of separable matrices to include the decomposition of matrices of rank < 2. This generalization leads to a straightforward polynomial time algorithm which computes 112
PAGE 119
113 this decomposition, i. e. the representation of a matrix of rank < 2 as the maximum of at most two outer products. 3) We presented a heuristic algorithm which computes the weak decomposition of matrices of arbitrary rank into outer products of row and column vectors. In most practical cases Â— i.e. whenever this decomposition is suited for the purpose of template decomposition Â— the number of these outer products is equal to the rank of the original matrix. 4) We derived the Minimax Eigenvector Decomposition (MED) within minmax algebra. In spite of its similarities to the Singular Value Decomposition, we showed that Â— unlike the SVD Â— the MED does not yield the best approximation of a matrix with respect to any specific norm. 5) Generalizing the problem of morphological template decomposition, we formulated the problem of morphological template approximation. We presented an integer programming approach to the problem of morphological template approximation. Clearly, this approach always guarantees the existence of a solution and can also be used for solving template decomposition problems. Moreover, slight variations in the original template will not lead to large variations in the objective function, which is the difference between the original matrix and the approximation with respect to a certain norm. 6) Justifying our integer programming approaches to the problem of morphological template approximation, we proved that this problem is NPcomplete. In fact, we showed that the following subproblem 11 is already NPcomplete: Given a matrix A and an integer k, is there a collection of column vectors u' and row vectors v', where
PAGE 120
114 / = l,...,k, such that A = \/ (u'm v') ? Further research is needed in order to answer the following questions: 1) What is the best strategy to solve the integer programming problems which we derived from the problem of morphological template approximation? 2) We presented a heuristic algorithm which can be used in order to compute the weak decomposition of a rectangular template into a minimal number of outer product expansions of column and row templates. Can this algorithm be generalized to include the decomposition of invariant templates into invariant templates having various supports? Can we formulate a similar heuristic algorithm which computes the decomposition of a rectangular template into 3x3 templates? 3) Can we prove the NPcompleteness of subproblems of the morphological template decomposition problems other than 11 ? For example, is the problem of strongly decomposing a rectangular template into 3x3 templates NPcomplete?
PAGE 121
REFERENCES [I] H.C. Andrews and B.R. Hunt. Digital Image Restoration. Prentice Hall, Reading, MA, 1977. [2] H.C. Andrews and C.L. Patterson. Outer product expansions and their uses in digital image processing. American Math. Monthly, 82(1): 113, January 1974. [3] R.C. Backhouse and B. Carr6. Regular algebra applied to pathfinding problems. J. Inst. Math. AppL, 15:161186, 1975. [4] D.H. Ballard and CM. Brown. Computer Vision. PrenticeHall, Englewood Cliffs, NJ, 1982. [5] C. Benzaken. Structures alg6bra des cheminements. In G. Biorci, editor, Network and Switching Theory, pages 4057. Academic Press, New York, 1968. [6] G. Birkhoff. Lattice Theory. American Mathematical Society, Providence, RI, 1984. [7] G. Birkhoff and J. Lipson. Heterogeneous algebras. J. Combinatorial Theory, 8:115133, 1970. [8] B. Carre. An algebra for network routing problems. J. Inst. Math. AppL, 7:273294, 1971. [9] J.E. Cohen. Subadditivity, generalized products of random matrices and operations research. SIAM Review, pages 6986, March 1988. [10] J.W. Cooley, P. A. Lewis, and P.D. Welch. The fast Fourier transform and its applications. IEEE Trans. Educ, E12(l):2734, 1969. [II] F. J. Crosby. Maxpolynomials and Morphological Template Decomposition. PhD thesis. University of Florida, Gainesville, FL, 1995. [12] H. Crowder, E.L. Johnson, and M.W. Padberg. Solving largescale zeroone linear programming problems. Operations Research, 31:803834, 1982. [13] R. CuninghameGreen. Process synchronisation in steelworks Â— a problem of feasibility. In Banbury and Maitland, editors, Procceedings of the 2nd International 115
PAGE 122
116 Conference on Operations Research, pages 323328, London, I960. English University Press. [14] R. CuninghameGreen. Minimax Algebra: Lecture Notes in Economics and Mathematical Systems 166. SpringerVerlag, New York, 1979. [15] J.L. Davidson. Lattice Structures in the Image Algebra and Applications to Image Processing. PhD thesis, University of Florida, Gainesville, FL, 1989. [16] J.L. Davidson. Minimax techniques for nonlinear image processing transforms. In Technical Symposium on Optics, ElectroOptics, and Sensors, volume 1098 of Proceedings of SPIE, Orlando, FL, March 1989. [17] J.L. Davidson. Foundations and applications of lattice transforms in image processing. In P. Hawkes, editor. Advances in Electronics and Electron Physics, pages 63130, New York, 1992. Academic Press. [18] J.L. Davidson. Nonlinear matrix decompositions and an application to parallel processing. Journal of Mathematical Imaging and Vision, l(2):2836, 1992. [19] J.L. Davidson. Classification of lattice transformations in image processing. Computer Vision, Graphics, and Image Processing: Image Uruler standing, 57(3):283306, May 1993. [20] P.D. Gader. Image Algebra Techniques for Parallel Computation of Discrete Fourier Transforms and General Linear Transforms. PhD thesis. University of Florida, Gainesville, FL, 1986. [21] P.D. Gader. Necessary and sufficient conditions for the existence of local matrix decompositions. SIAM Journal on Matrix Analysis and Applications, pages 305313, July 1988. [22] P.D. Gader. Separable decompositions and approximations for grayscale morphological templates. CGVIP, 53:288296, July 1991. [23] P.D. Gader and S. Takriti. Decomposition techniques for grayscale morphological templates. SPIE, 1350:431^2, July 1990. [24] P.D. Gader and S. Takriti. Local decomposition of grayscale morphological templates. Journal of Mathematical Imaging and Vision, 2:3950, 1992. [25] M.R. Garey and D.S. Johnson. Computers and Intractibility: A Guide to the Theory of NPCompleteness. W.H. Freeman and Company, New York, 1979. [26] A.M. Geoffrion and R.E. Marsten. Integer programming algorithms: A framework and a state of the art survey. Management Science, 18:465491, 1972.
PAGE 123
117 [27] B. Giffler. Mathematical solution of production planning and scheduling problems. Tech. rep., IBM ASDD, 1960. [28] G.H. Golub and C.F. van Loan. Matrix Computations. John Hopkins University Press, Baltimore, MD, 1983. [29] R.C. Gonzalez and P. Wintz. Digital Image Processing. AddisonWesley, Reading, MA, second edition, 1987. [30] R.C. Gonzalez and R.E. Woods. Digital Image Processing. AddisonWesley, Reading, MA, 1992. [31] H. Hadwiger. Vorlesungen Uber Inhalt, Oberflceche und Isoperimetrie. SpringerVerlag, Berlin, 1957. [32] H. J. A. M. Heijmans. Morphological Image Operators. Academic Press, Boston, MA, 1994. [33] H. J. A. M. Heijmans. Mathematical morphology: A modem approach in image processing based on algebra and geometry. SI AM Review, 37(1): 136, 1995. [34] N. Karmarkar and K. Ramakrishnan. Computational results of an interior point algorithm for large scale linear programming. Mathematical Programming, 52:555585, 1991. [35] A. K. Katsaggelos. Digital Image Restoration. Springer Verlag, Berlin, 1991. [36] S.Y. Lee and J.K. Aggarwal. Parallel 2D convolution on a mesh connected array processor. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI9(4):590594, July 1987. [37] D. Li. Morphological template decomposition with maxpolynomials. Journal of Mathematical Imaging and Vision, 1(3):215221, September 1992. [38] D. Li and G. X. Ritter. Decomposition of separable and symmetric convex templates. In Image Algebra and Morphological Image Processing, volume 1350 of Proceedings of SPIE, pages 408^18, San Diego, CA, July 1990. [39] Z.Z. Manseur and D.C. Wilson. Decomposition methods for convolution operators. Computer Vision, Graphics, and Image Processing, 53(5):428^34, 1991. [40] P. Maragos and R.W. Schafer. Morphological skeleton representation and coding of binary images. IEEE Trans. Acoustics, Speech,and Signal Proc, ASSP34(5):12281244, October 1986. [41] G. Matheron. Random Sets and Integral Geometry. Wiley, New York, 1975.
PAGE 124
118 [42] H. Minkowski. Volumen und oberflache. Mathematische Annalen, 57:447^95, 1903. [43] H. Minkowski. Gesammelte Abhandlungen. Teubner Verlag, LeipzigBerlin, 1911. [44] D.P. O'Leary. Some algorithms for approximating convolutions. Computer Vision Graphics and Image Processing, 41(3):333345, 1988. [45] P.M. Pardalos and Y. Li. Integer programming. In C.R. Rao, editor. Handbook of Statistics, volume 9, pages 279302. Elsevier, 1993. [46] H. Park and R. T. Chin. Optimal decomposition of convex morphological structuring elements for 4connected parallel array processors. IEEE Trans, on Pattern Anal, and Mach. Intell, 16(3):304313, 1994. [47] G.X. Ritter. Recent developments in image algebra. In P. Hawkes, editor. Advances in Electronics and Electron Physics, volume 80, pages 243308. Academic Press, New York, NY, 1991. [48] G.X. Ritter. Image Algebra, unpublished manuscript, available via anonymous ftp from ftp.cis.ufl.edu in /pub/ia/documents, 1995. [49] G.X. Ritter. Handbook of Computer Vision Algorithms in Image Algebra. CRC Press, Boca Raton, FL, 1996. [50] G.X. Ritter and P.D. Gader. Image algebra: Techniques for parallel image processing. Journal of Parallel and Distributed Computing, 4(5):744, April 1987. [5 1 ] G.X. Ritter and P. Sussner. The minimax eigenvector transform. In Image Algebra and Morphological Image Processing III, Proceedings of SPIE, San Diego, CA, July 1992. [52] G.X. Ritter, J.N. Wilson, and J.L. Davidson. Image algebra: An overview. Computer Vision, Graphics, and Image Processing, 49(3):297331, March 1990. [53] L. Schrage. Linear, Integer and Quadratic Programming with UNDO. The Scientific Press, Palo Alto, CA, 1983. [54] J. Serra. Image Analysis and Mathematical Morphology. Academic Press, London, 1982. [55] J. Serra. Image Analysis and Mathematical Morphology, Volume 2: Theoretical Advances. Academic Press, New York, 1988. [56] S.R. Sternberg. Biomedical image processing. Computer, 16(1), January 1983.
PAGE 125
119 [57] P. Sussner, P. M. Pardalos, and G. X. Ritter. Global optimization problems in computer vision. In C. Floudas and P. Pardalos, editors. State of the Art in Global Optimization, pages 457474. Kluwer Academic Publishers, 1996. [58] S. Takriti and P.D. Gader. Local decomposition of grayscale morphological templates. Journal of Mathematical Imaging and Vision, 2:3950, 1992. [59] University of Florida Center for Computer Vision and Visualization, Gainesville. Handbook of Image Algebra, 1993. [60] J. Xu. Decomposition of convex polygonal morphological structuring elements into neighborhood subsets. IEEE Transactions on Pattern Analysis and Machine Intelligence, 13(2), February 1991. [61] H. Zhu and G. X. Ritter. The pproduct and its applications in signal processing. SIAM Journal on Matrix Analysis and Applications, 16(2):579601, April 1995. [62] X. Zhuang and R. M. Haralick. Morphological structuring element decomposition. Computer Vision, Graphics, and Image Processing, 35:370382, 1986.
PAGE 126
BIOGRAPHICAL SKETCH Peter Sussner was bom in Nuremberg, Germany. He studied mathematics and computer science at the University of Erlangen. He graduated with a M. S. in mathematics in 1990. One year later, he continued his studies at the University of Florida, where he became a research assistant at the Center for Computer Vision and Visualization. His research interests include applied mathematics, artificial neural networks, image processing, and computer vision. 120
PAGE 127
I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in sco^ and quality, as a dissertation for the degree of Doctor of Philosophy. Qgward Ritter, Chairman rofessor of Mathematics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. William Hager Professor of Mathematics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Panos Pardalos Associate Professor of Industrial and Systems Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Murali Rao Professor of Mathematics I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. David Wilson Professor of Mathematics
PAGE 128
I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. This dissertation was submitted to the Graduate Faculty of the Department of Mathematics in the College of Liberal Arts and Sciences and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. August 1996 ^oseplr^ilson Assistant Professor of Computer and Information Science and Engineering Dean, Graduate School

