The Piecewise Linear-QuadraticModel for Convex BivariateFunctionsbyScott Ronald FazackerleyA THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENTS FOR THE DEGREE OFBachelor of Science HonoursinThe I. K. Barber School of Arts & Sciences(Computer Science)The University Of British Columbia(Kelowna, Canada)April, 2008a169 Scott Ronald Fazackerley 2008AbstractA Piecewise Linear-Quadratic (PLQ) function is used to describe data that can be repre-sented by continuous functions with a piecewise linear domain for which the function is eitherlinear or quadratic on, each piece of its domain [11]. While extensive work has been donewith PLQ models for convex univariate functions, my work investigates the development ofa two dimensional model that allows the implementation of algorithms for computing funda-mental convex transforms. PLQ functions have been described in the literature, the efficientimplementation of the algorithms requires the careful selection of a data structure.Initial investigation examined using a Voronoi diagram to represent the projection of thebivariate function into a822, to take advantage of efficient algorithms for the point locationproblem [3]. While suitable for representing a single PLQ, with operations for building a zeroorder model from data and evaluating the function over an irregular grid, we are uncertainif it can be extended to compute other fundamental convex transforms efficiently.In examining the Voronoi model, we were able to extend the base concept to represent thebivariate PLQ using two different representations: a tessellation-based model and a linearinequality-based model. The tessellation model is restricted to representing a PLQ with abounded domain. The model represents data through triangular faces in a tessellation in a822where each face is defined by its vertices and the associated function value. This model allowsfor the efficient evaluation over an irregular grid, the addition of two PLQ functions withbounded domains and multiplication by a scalar, but does not allow for the representationof an unbounded domain. To allow for an unbounded domain, a dual model was developedusing linear inequalities to represent each face in a822. Numerical results are presented foreach model and computational complexity of model components are discussed.iiAcknowledgmentsI would like to thank Dr. Yves Lucet for providing me with the opportunity to do an Honoursunder his supervision. His feedback and insights have been invaluable in bringing my thesisto its current form. This has been an exceptional learning experience in areas of computingand mathematics which I would not have otherwise had the opportunity to explore.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Model Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.1 Tessellation-based Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.2 Inequality-based Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Algorithms for the Bounded Tessellation-Based Bivariate PLQ Model . 103.1 Build Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103.2 Evaluation Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113.3 Addition Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133.4 Scalar Multiplication Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 144 Algorithms for the Unbounded Inequality-Based Bivariate PLQ Model 154.1 Algorithm for Converting a Tessellation-Based Model to an Inequality-BasedModel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.2 Evaluation Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 204.3 Addition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224.3.1 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 224.4 Scalar Multiplication Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 245 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . 265.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26iv6 APPENDIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 296.1 Building a Tessellation-Based PLQ from Pointwise Data . . . . . . . . . . . 296.1.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 296.1.2 Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 316.2 Evaluating a Bounded PLQ . . . . . . . . . . . . . . . . . . . . . . . . . . . 316.2.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 316.2.2 Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 326.3 Adding Two Bounded PLQ’s . . . . . . . . . . . . . . . . . . . . . . . . . . 326.3.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 326.3.2 Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 336.4 Scalar Multiplication of a tessellation-based PLQ . . . . . . . . . . . . . . . 356.4.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 356.4.2 Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 356.5 Converting a tessellation-based PLQ to an inequality-based PLQ . . . . . . 376.5.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 376.5.2 Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 386.6 Evaluating an inequality-based PLQ . . . . . . . . . . . . . . . . . . . . . . 396.6.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 396.6.2 Comments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 396.7 Additional Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41vList of Figures1.1.1 Related Voronoi and Delaunay Diagrams . . . . . . . . . . . . . . . . . . . . 31.2.1 Model Representation over a Convex Set . . . . . . . . . . . . . . . . . . . . 53.2.1 Empirical results from the tsearch algorithm in MATLAB . . . . . . . . . . . 136.1.1 Build Method Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 306.1.2 Evaluation of a tessellation-based PLQ2 . . . . . . . . . . . . . . . . . . . . 316.3.1 The Addition of Two PLQ Models . . . . . . . . . . . . . . . . . . . . . . . 346.4.1 PLQ Multiplication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 366.5.1 Run Time Complexity for Model Conversion . . . . . . . . . . . . . . . . . . 376.6.1 IPLQ Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40viList of Tables4.1.1 Sign value encoding for inequalities for the bounded polytope . . . . . . . . . 164.1.2 Sign value encoding for inequalities for bounded and unbounded regions . . . 175.1.1 Complexity Results for Fundamental Convex Transforms . . . . . . . . . . . 266.1.1 Complexity Results of the Build Method for a Strictly Convex Function . . . 31viiChapter 1IntroductionExtensive work has been completed in regards to a model for a univariate Piecewise Linear-Quadratic (PLQ) function. A univariate PLQ is used to describe data that can be representedby continuous functions with a piecewise linear domain for which the function is either lin-ear or quadratic on each piece of its domain [11]. Univariate PLQ functions appear in thefields of economic modeling, control theory and scientific computing, specifically computa-tional convex analysis [10, 11, 13]. While work has been done with PLQ models for convexunivariate functions, a bivariate PLQ model has not yet been presented that supports ef-ficient implementation of fundamental convex transforms. PLQ functions are closed understandard convex operations, and can be manipulated with linear-time algorithms [11]. Effi-cient manipulation is defined as having a storage and run-time complexity that is at worstpolynomial.My work investigates the development of a two dimensional model that allows the imple-mentation of algorithms for computing fundamental convex transforms. While PLQ func-tions have been described in the literature, the efficient implementation of the algorithmsrequires the careful selection of a data structure. The requirements for the model are pre-sented in Section 1.2. Three models are examined and discussed in Section 2.The models presented are zero-order models: only point data can be accepted as input tothe model. The algorithms presented in Section 3 and Section 4 compute function values overthe domain of the input data. Convexity is enforced over the domain of input values: pointsmaking the function non-convex are removed such that the function remains continuousand convex over the domain. Additionally, algorithms for converting from a bounded to anunbounded model, evaluating, adding and (scalar) multiplying are discussed. Complexityresults are also presented, proven and validated with empirical data.1.1 PreliminariesIn this section, we recall some basic notions from computational geometry and convex anal-ysis.Definition 1.1.1. Convex Set [6] A set C is convex, if the line segment between any twopoints in C lies in C. That is, for any x1,x2 ∈ C and any θ with 0 ≤ θ ≤ 1, we haveθx1 + (1−θ)x2 ∈ C. (1.1.1)Definition 1.1.2. Convex Hull [6] The convex hull of a set C, denoted conv C is the setof all convex combinations of points in C1conv C = {θ1x1 + ... + θkxk|xi ∈ C,θi ≥ 0,i = 1,...,k,θ1 + ... + θk = 1}. (1.1.2)The convex hull conv C is always convex.Definition 1.1.3. Convex Function [6] A function f : Rn → R ∪{+∞} is convex if domf = {x|f(x) < ∞} is a convex set and if for all x, y ∈ domf, and θ with 0 ≤ θ ≤ 1f(θx + (1−θ)y) ≤ θf(x) + (1−θ)f(y). (1.1.3)A function f is strictly convex if strict inequality holds in Equation (1.1.3) for every x negationslash= yand 0 < θ < 1.Definition 1.1.4. Triangulation [8] A triangulation is a planar subdivision whose boundedfaces are triangles. We note P := {p1,p2,...,pn} is the set of points in the plane forming thetriangulation.Definition 1.1.5. Maximal Planar Subdivision [8] A maximal planar subdivision is a subdi-vision S such that no edge connecting two vertices can be added to S without destroying itsplanarity.Definition 1.1.6. Euclidean Distance [8, pp. 148–149] The Euclidean distance between twopoints p and q is denoted as dist(p,q). Thus, in the plane we havedist(p,q) =radicalBig(px −qx)2 + (py −qy)2. (1.1.4)Definition 1.1.7. Voronoi Diagram [8, pp. 148–149] Let P := {p1,p2,...,pn} be a set of ndistinct points in the plane; these points are the sites. The Voronoi diagram of P is definedas the subdivision of the plane into n cells, one for each site in P, with the property that apoint q lies in the cell corresponding to a site pi if and only if dist(q,pi) < dist(q,pj) for allpj ∈ P with i negationslash= j. The Voronoi diagram of P is denoted by Vor(P). The region of a sitep is called the Voronoi cell of p and is denoted by V(p). A Voronoi diagram is pictured inFigure 1.1(a).Definition 1.1.8. Delaunay Triangulation [8, pp. 190–191] Let the graph G be the dualgraph of a Voronoi diagram constructed as follow: The graph G has a node for every Voronoicell and it has an arc between two nodes if the corresponding cells share an edge. This meansthat G has an arc for every edge of Vor(P). There is a one-to-one correspondence betweenthe bounded faces of G and the sites of Vor(P). Consider the straight-line embedding of G,where the node corresponding to the Voronoi cell V(p) is the point p, and the arc connectingthe nodes of V(p) and V(q) is the segment pq. This embedding is called the Delaunay graphof P and is denoted by DG(P). A Delaunay triangulation is any triangulation obtained byadding edges to the Delaunay graph. The Delaunay triangulation of P is unique if and onlyif DG(P) is a triangulation. Thus a Delaunay triangulation can be characterized by:2(a) Voronoi Diagram (b) Delaunay TriangulationFigure 1.1.1: Related Voronoi and Delaunay DiagramsLet P be a set of n points in the plane.i Three points pi,pj,pk ∈ P are vertices of the same face of the Delaunaygraph of P if and only if the circle through pi,pj,pk contains no pointsof P in its interior.ii Two points pi,pj ∈ P form an edge of the Delaunay graph of P if andonly if there is a closed disc C that contains pi and pj on its boundaryand does not contain any other point of P.A Delaunay triangulation is pictured in Figure 1.1(b). This is the dual structure to theVoronoi diagram pictured in Figure 1.1(a).Definition 1.1.9. The centroid of a finite set of points p1,p2,..pn is the arithmetic mean1nnsummationdisplayi=0p1. (1.1.5)It lies inside the convex hull [12, p. 106].Definition 1.1.10. Piecewise Linear-Quadratic (PLQ) Function [11] A function is a PLQfunction if it can be represented by continuous functions with a piecewise linear domain forwhich the function is either linear or quadratic on each piece of its domain.1.2 Problem DefinitionWe desire to produce a model for a bivariate PLQ function that will be both efficient in spaceand run-time complexity for building and representing the PLQ as a numerical model, as3well as having an efficient run-time complexity under fundamental convex transforms. Theinitial goal is to represent zero order data, with the user providing input data in terms of aseries of points as (xi,yi,zi) in a823. The model is required to remove any co-planar pointsinside a face on the convex hull. This results in storage complexity of the model being lessthan or equal to that of the size of the initial data set. This will minimize the amount ofstorage required to represent the model without loss of generality.Figure 1.1(a) represents 401 input points over a regular grid. When represented with thePLQ model, as there are numerous co-planar points, the results are simplified such that themodel is defined by only 9 vertices as in Figure 1.1(b). With a strictly convex input set as inFigure 1.1(c), which represents 401 input points over a regular grid, no simplification occursin the representation of the model. Figure 1.1(d) contains 401 vertices due to the fact thatthe function is strictly convex and no simplification can occur. With the zero-order model,limitations in terms of complexity will be encountered with data that represents a strictlyconvex function as no more than three points are co-planar on the convex hull. Thus nosimplification will occur, presenting a upper bound for storage complexity.In addition to supporting algorithms to allow for the efficient transformation of the modelfrom one form to another, the model must support the following operations in an efficientmanner:a136 building a model from zero-order data,a136 evaluating the model over an irregular grid of points,a136 multiplying by a scalar value, anda136 adding of two PLQ functions.4(a) Convex input points for f(x,y) = |x|+|y| (b) Zero-Order model of f(x,y) = |x|+|y|(c) Strictly Convex input points for f(x,) = x2+y2 (d) Zero-Order Model of f(x,y) = x2 + y2Figure 1.2.1: Model Representation over a Convex Set5Chapter 2Model Selection2.1 Tessellation-based ModelIn examining the problem, it was felt that a PLQ function could be encoded using a Delaunaytriangulation, where each data point in the zero-order model would represent a vertex. TheDelaunay triangulation is a dual structure of the Voronoi diagram, which has previouslybeen used to encode a convex hull as a finite set of intersecting half-spaces [3], where eachface is defined by a Voronoi center V.The Delaunay triangulation was chosen over the Voronoi diagram as we could explicitlydefine the domain over which a function is defined. A user would define the vertex pointvalue in the model and the algorithm would compute and project triangulations into a822.Encoded with each triangulation would be the value of the function as a linear-quadraticequation. Thus a PLQ function in a823 is encoded as a matrix containing triangulations as Tin a822 associated with the corresponding linear-quadratic equation asf(Xi) =parenleftbigg ai,1 ai,2ai,3 ai,4parenrightbiggX2i +parenleftbigg bi,1bi,2parenrightbiggXi + ci. (2.1.1)where Xi is the point (xi,yi) at which the linear-quadratic function is to be evaluated.Equation (2.1.1) encodes the function value of each part of the PLQ to the boundaries ofthe polytope. To encode the information as a one-dimensional array we rewrite Equation(2.1.1) as parenleftbigai1,1 ai1,2 ai2,1 ai2,2 bi1,1 bi1,2 ci parenrightbig. (2.1.2)By definition, the points P := {p1,p2,...,pn} that form the set of points in the plane,have a one-to-one mapping onto T , but now we can explicitly define the boundaries overwhich Equation (2.1.1) applies. What we need is not unlike the process of building a terrainmap, but with the constraints applied by the conditions of a bivariate PLQ.The model that we propose, is a function proj : A ⊂ a823 → a822, that assigns a linear-quadratic function value, f(p) to every point within a planar subdivision in the domain A.As with a terrain map, the value of the bivariate PLQ is not necessary known at every pointin the zero-order model, but defined at certain points such that we only know the value ofthe function proj at a finite set P ⊂ A [8].The initial challenge is how to determine a triangulation of P, such that each pi ∈ P is avertex for each planar sub-division and that each planar sub-division is a triangle. Followingthat, we then need to project the triangulation from a823 → a822. Additionally, we needto compute the function coefficient values for Equation (2.1.1) for each planar subdivision.6Initial investigation was completed with the Scilab mathematical package [2]. Unfortunately,Scilab does not have the necessary functionality to support the operations we were lookingfor without a significant amount of time invested into writing the basic algorithms. Themathematical package MATLAB [1] was found to offer the necessary functionality to supportthe desired operations. The nature of the specific algorithms will be expanded in the followingsections.Using MATLAB as the computational platform, we were able to develop a data structurethat met the requirements to encode a bivariate PLQ. Each planar subdivision is defined by atessellation which is composed of two separate parts: one being the set of points (xi,yi) ∈ P,which is now a maximal planar subdivision, and the other being a triangulation of P, whichis defined as followsDefinition 2.1.1. Triangulation of P [8, p. 185] A triangulation T of P is a maximalplanar subdivision whose vertex set is PTo build the bivariate PLQ from zero-order data, a convex hull algorithm was run on thedata set ina823 . The resulting convex hull was then projected intoa822, providing the requiredtriangulations. Additionally, the coefficient values for Equation (2.1.1) were computed fromvertex points in the convex hull, forming the necessary components to represent the bivariatePLQ. The data structure has two distinct components: a matrix of (xi,yi) ∈ P where Pcontains the location in a822 of each vertex. A second matrix encodes triangulation ti ∈ D asa set of indexes into P along with the corresponding linear-quadratic as in Equation (2.1.2).Let D be a triangulation T that satisfies the properties of a Delaunay tessellation. Themodel can be writtenP =x1 y1x2 y2· ·· ·· ·xnv ynv(2.1.3)D =i1,1 i1,2 i1,3 a1,1 a1,2 a1,3 a1,4 b1,1 b1,2 c1i2,1 i2,2 i2,3 a2,1 a2,2 a2,3 a2,4 b2,1 b2,2 c2· · · · · · · · · ·· · · · · · · · · ·int,1 int,2 int,3 ant,1 ant,2 ant,3 ant,4 bnt,1 bnt,2 cnt(2.1.4)where ij,k,(k ∈ {1,2,3} and j = 1..nt) are the indexes of the associated points in P with nvand nt being the number of vertexes and triangulations respectively.This model supports the desired convex transformation operations efficiently. The de-scription and analysis of each component is discussed in Chapter 3. The storage complexityof the model is O(n) in terms of the number of triangulations as there is one entry for eachtriangulations and at most 3 vertexes for each triangulation. The only shortcoming to rep-resenting a bivariate PLQ with this data structure is that it cannot represent an unboundedbivariate PLQ due to the fact that the union of the bounded faces of a triangulation D is7the convex hull of P [8]; hence, it cannot encode an unbounded region and is limited torepresenting bounded domains.2.2 Inequality-based ModelIn order to support the inclusion of unbounded domains into the model, a dual model wasinvestigated using inequalities instead of vertices in P to define each triangulation ti ∈ T .This is accomplished through a series of operations that convert a tessellation-basedPLQ to an inequality-based PLQ. Each ti ∈ T is now represented by several half-spacesin a822 with the associated linear-quadratic equation. The matrix P is replaced by a listof inequalities I, where (mi,si,bi) ∈ I are the set of half-spaces that define T . Eachcomponent of (mi,si,bi) ∈ I encodes the slope m, the sign of the inequality s and theintercept b respectively. Let I be a set of inequalities, T be a set of triangulations, C be aset of centroids and N be a set of neighbouring faces. The model becomesI =m1 s1 b1m2 s2 b2· · ·· · ·· · ·mns sns bns(2.2.1)T =i1,1 i1,2 i1,3 a1,1 a1,2 a1,3 a1,4 b1,1 b1,2 c1i2,1 i2,2 i2,3 a2,1 a2,2 a2,3 a2,4 b2,1 b2,2 c2· · · · · · · · · ·· · · · · · · · · ·· · · · · · · · · ·int,1 int,2 int,3 ant,1 ant,2 ant,3 ant,4 bnt,1 bnt,2 cnt(2.2.2)C =x1 y1x2 y2· ·· ·· ·xnt ynt(2.2.3)N =i1,1 · i1,nn1i2,1 · i2,nn2· · ·· · ·· · ·int,1 · int,nnt(2.2.4)where ij,k,k ∈ {1,2,3} in the second matrix are the indexes of the associated inequalities inI and ns, nt, nnt are the number of half-spaces, the number of triangulations and the number8of neighbours for a specific triangulation respectively. The array of lists N maintains theindexes of each neighbouring triangulation is T to a triangulation ti.The inequality model introduces two additional components that are used to support theefficient transformations as a third and fourth element. The third element is a matrix ofcentroids. Each row stores the centroid for the triangulation defined by the given row. Thefourth element is a cell list, and stores a list of indexes to the neighbours of the triangulationdefined by a given row. A cell list is a data structure specific to MATLAB which allowseach row to have a different number of elements. This is critical and is done to minimize thestorage complexity as each triangulation can have a different number of neighbour. This listmust also be able to grow in case new neighbouring triangulations are added.This model supports the desired convex transformation operations in terms of efficientevaluation, addition and multiplication, but must be constructed from a bounded bivariatePLQ model. The description and analysis of each component is discussed in Chapter 4.9Chapter 3Algorithms for the BoundedTessellation-Based Bivariate PLQModel3.1 Build AlgorithmIn constructing a bivariate PLQ, both the tessellation and the inequality-based models relyon the same initial build function. The build function, working only with zero-order data,receives as input an irregular grid of points X in a822 and a vector of corresponding functionvalues, f(X). The build algorithm isAlgorithm 1: PLQ2 build AlgorithmInput: X,fOutput: a tessellation-based bounded PLQ2begininitialization;Compute the convex hull, H of the data set in a823;Project the vertices P that form H into a822;Compute the Delaunay triangulation D ∀p ∈ P;forall ti ∈ D docompute the coefficients for linear-quadratic function;endReturn PLQ;endThe convex hull algorithm used in this implementation, is based on the Quickhull algo-rithm [4]. The algorithm will remove extra points that are located on the interior of eachtriangulation [8, 4].Let us recall some factsa136 The number of points processed by the convex hull algorithm is proportional to thenumber of vertices in the output and runs in O(nlogh), where h is the output size [7].a136 A convex polytope P with n vertices has at most 2n−4 faces [8, p. 237].a136 The implementation of Quickhull in MATLAB produces a convex hull with triangularfaces with three vertices.10a136 The expected number of triangles created by the Delaunay Triangulation algorithm isas most 9∗n + 1 [8, p. 197].Proposition 3.1.1. The PLQ2 build algorithm runs in O(nlogh).Proof. Let P be a convex polytope produced by the algorithm in O(nlogh) which has atmost 2h−4 facets. The projection of P from a823 →a822 requires each face to be visited onlyonce and for each face, only three points are projected. Thus, let P be the set of points ina822 and the total number of operations to project into a822 is 3∗(2n−4) = 6n−12 = O(n).By Definition 1.1.8, the Delaunay triangulation D of a set P of 6n − 12 points can becomputed in O(nlogn) and produces at most 9 ∗ (6n − 12) − 1 = 54n − 107 faces. For asingle face, the linear-quadratic equation coefficients is computed in O(1). It then followsthat the coefficients for all triangulations in D can be computed in O(54n−107) = O(n). Sothe total number of operations required for the PLQ2 build algorithms is O(nlogh)+O(n)+O(nlogn) + O(n) = O(nlogh).3.2 Evaluation AlgorithmThe evaluation of the bounded bivariate PLQ function over a set of points (xi,yi) ∈ X isachieved exploiting the fact that the model is based on a tessellation. MATLAB includesa function called tsearch in the convex hull package that allows a point or set of points tobe evaluated over a Delaunay triangulation D and will return the index of the face that thepoint is included in. By definition D encodes a convex hull which is a continuous piecewisefunction. If a point (x,y) falls on the boundary of a triangulation, it can be evaluated aspart of either triangulation. If a point that is being evaluated is not in the domain of themodel, the value returned is +∞.Algorithm 2: PLQ2 eval AlgorithmInput: PLQ2,XOutput: a vector of values Zbegininitialization;forall (xi,yi) ∈ D doCompute facet indexes j ∈ I from D;endforall i doevaluate the linear-quadratic function at (xi,yi) associated with the facet atindex ji;insert results of evaluation at Z(i);endReturn Z;endThe computational complexity of the tsearch algorithm as published in MATLAB wasnot available through the literature. Thus, an empirical approach was taken to evaluate11the run-time complexity of the bounded evaluation algorithm. A series of experiments wererun with two different scenarios: the first being evaluating a fixed number of points over anincreasing number of facets in the model. The second being the evaluation of an increasingnumber of points over a fixed number of facets in the model.Conjecture 3.2.1. Tsearch runs in O(nk), where n in the number of facets in D and k isthe number of points being evaluated.In Figure 3.1(a), 100 points were evaluated in a PLQ model that had a maximal numberof n facets (strictly convex data) with n growing to 80,000. It was observed that the timefor evaluating a fixed number of points increased linearly with the number of facets in themodel. In Figure 3.1(b), the number of input points to be evaluated was increased up to250,000 points over a model with 80,000 facets in D. The computational time was foundto increase linearly as the number of points being evaluated. Thus the tsearch algorithmappears to run in O(nk).Returning to the PLQ2 eval algorithm, we continue with analyzing the remainder of thealgorithm.Proposition 3.2.1. Assuming that tsearch runs in O(nk), the PLQ2 eval algorithm runs inO(nk).Proof. Let X be a set of input points of size k over D with n facets. From Conjecture 3.2.1,the runtime complexity will be O(nk) producing an index list I of size k as per the definitionof the function in MATLAB. To evaluate the function value at each point in X, the linear-quadratic function to be evaluated can be accessed and evaluated in O(1) + 7 = O(1). Itthen follows that accessing k linear-quadratic functions can be completed in O(k) and theruntime complexity of the algorithm is O(nk) + O(k) = O(nk).Note: A better implementation of tsearch may improve its complexity to O(nlogk)therefore improving the complexity of PLQ2 eval to O(nlogk + k).12(a) Timing results for evaluating a fixed number ofpoints in the tessellation over an increasing numberof facets(b) Timing results for evaluating a increasing num-ber of input points over 80,000 facetsFigure 3.2.1: Empirical results from the tsearch algorithm in MATLAB3.3 Addition AlgorithmThe PLQ2 add algorithm takes as input, two bounded bivariate models and returns a bi-variate PLQ that is the addition of the two input models. The addition algorithm is a unionof the vertex points P in D under evaluation of both input models. The algorithm is asfollows.Algorithm 3: PLQ2 add AlgorithmInput: PLQa,PLQbOutput: a PLQ that is a summation of the two input modelsbegininitialization;for each vertex point pi ∈ Pa in PLQ2a doevaluate the function value of pi under PLQb to form a vector Za;endfor each vertex point pj ∈ Pb in PLQ2b doevaluate the function value of pj under PLQa to form a vector Zb;endConcatenate Pa and Pb to form a new point list Pc;Concatenate Za and Zb to form a new point list Zc;Evaluate Pc and Zc under the PLQ2 build to build PLQ2c;Return PLQ2c;endProposition 3.3.1. The PLQ2 add algorithm runs in O(n2).13Proof. Let n be the number of triangulations in PLQa and m be the number of triangulationsin PLQb. Let j be the number of vertex points in PLQa (at most j = 3n) and k be thenumber of vertex points in PLQb (at most k = 3m). It then follows that the evaluationof PLQa in PLQb is O(3nm) and the evaluation of PLQb in PLQa is O(3mn). Whenn = m, the total cost is O(3n2) + O(3n2) = O(n2) with j + k = 3m + 3n = 6n points.The computational complexity of building the resulting new PLQ from the data pointsby Proposition 3.1.1 is O(6nlog6n) = O(nlogn). It then follows that the computationalcomplexity for the addition algorithm is O(n2) + O(nlogn) = O(n2).3.4 Scalar Multiplication AlgorithmScalar multiplication for a PLQ function is a straightforward linear-time algorithm in whichthe linear-quadratic equation for each face in the PLQ model is multiplied by a scalar value.Algorithm 4: PLQ2 mult AlgorithmInput: PLQ2,cOutput: a PLQ2 that has been scalar multiplied by cbegininitialization;forall triangulations ti ∈ D in PLQ2 domultiply the linear-quadratic equation at index i by the constant c;endReturn PLQ2;endProposition 3.4.1. The PLQ2 mult algorithm runs in O(n).Proof. Let n be the number of triangulations in the PLQ and let x be a scalar value. Loop-ing on n, indexing the value of the linear-quadratic function for each triangulation in R2can be completed in O(1) by the fact that indexing into a matrix can be done in lineartime. Computing the new linear-quadratic equation from the scalar multiplication of kwith the linear-quadratic equation, Equation (2.1.2), is completed in O(1). Updating thelinear-quadratic equation at index i can be done in linear time. In then follows that eachtriangulation’s linear-quadratic equation must be accessed taking n operations. Hence thetotal number of operations is n∗(O(1) + O(1) + O(1)) = O(n).14Chapter 4Algorithms for the UnboundedInequality-Based Bivariate PLQModel4.1 Algorithm for Converting a Tessellation-BasedModel to an Inequality-Based ModelThe algorithm to convert a bounded tessellation-based bivariate PLQ to an unboundedbivariate inequality-based PLQ takes as an argument a bounded PLQ and returns and un-bounded PLQ as IPLQ2. The algorithm is composed of three distinct operations: first, thebounded tessellation model is converted to an inequality model. Second, the unboundedregion is constructed as a IPLQ2 model, and third, the two models are joined together.The first set of operations is computed by treating each vertex point of a specific tri-angulation as the intersection of two inequalities formed by the intersection point and theremaining two vertices. This method exploits the property of the Delaunay Triangulation inMATLAB, in such that each face is triangular and will only have three vertices. The slopem and intercept value b of each inequality are computed from the vertex information. Eachinequality ineq is encoded in the point-slope form of y = mx + b as (m sign b). The in-equality value sign must then be determined. By the properties of a closed convex polytope,any point that is on the boundary of the polytope is within the convex hull, thus as eachinequality ineq is uniquely formed by only two vertices. The third vertex that is not takingpart in the formation of the inequality ineq is used as a test point as it will, by definition,be in the half-space that forms part of the initial triangulation. Values of sign are encodedas in Table 4.1.1.Special consideration must be given to cases where the slope m of the line is +∞ or −∞as the inequality is in the form of y = mx + b. To capture and encode this information, aspecial case is allowed where m can take on the values of +∞ or −∞. In this case, the valueof b encodes the value of the x-intercept.Inequalities, as described in Equation (2.2.1) are stored in a matrix of inequalities andeach face then maintains a list of indices to the required inequality. As the model is beingconstructed, duplicate checking is conducted to ensure that existence of an equality in thematrix is unique.The unbounded region around the convex polytope that represents the triangulation ina822 is computed by touring around the exterior of the convex hull. Based on the properties of15Sign Value≥ encoded as 1= encoded as 0≤ encoded as -1Table 4.1.1: Sign value encoding for inequalities for the bounded polytopethe convex hull, for a convex hull tour, the first point past the current segment being tested,will be used as a test point to determine the unbounded regions. The computation of theconvex hull of the bounded model is performed using the convhull function in MATLAB [1, 4]and the vertex points from the Delaunay triangulation is a822. This returns the vertices thatform the convex hull in a strict ordering, such that indexing through the vertices will form atour around the convex hull returning to the original vertex. The assumption is made thatthe hull consists of at least three points. If the convex hull has only two points, then it is aline and is treated as a special case.Unlike the bounded region where a point on the boundary of a polytope is inclusive ofthe region, in the unbounded region a point on the boundary of the convex polytope maynot belong to the half-space defined by the inequalities defining the region. Thus < and >were added as special cases for unbounded conditions: if a point is sitting on the convex hullof the tessellation in a822, it is in the bounded model. As there could exist a discontinuitybetween the convex set and the unbounded region, a point on the hull cannot be equal tothe linear-quadratic equation for both spaces. Thus, the unbounded region as expressed bythe inequalities, must bound the convex hull with a strict inequality, hence the need for theintroduction of < and > as in table 4.1.2.In the calculation of the inequalities forming the unbounded region, the tour around thehull will form a series of unbounded faces consisting of two inequalities; The first, will be thestrict inequality that forms the boundary of the convex hull. The second point may eitherbe an inequality or strict inequality as it is only forming a partition between two unboundedfaces that evaluate to +∞. Each unbounded region will be formed by no more than twoseparate inequalities. Again, as with the bounded region, cases of infinite slope are encodedin the same fashion. The tessellation of the unbounded regions are encoded in the samefashion with the IPQ2 data structure.With conversion of the bounded model to an unbounded model, being represented bytwo separate parts, the unbounded and bounded inequality models are merged into a singledata structure. This can be accomplished by a concatenation of the two models, but asat least half of the inequalities that form the unbounded region take part in the formationof the bounded inequality model, duplication of inequalities is inevitable. Thus, duplicateinequalities are removed and existence of inequalities in one list need to be verified in theother.MATLAB does not contain a native binary search method, which would be a suitablealgorithm for checking the existence of an inequality in a matrix of inequalities. In theMATLAB file exchange, a binary search algorithm (that given a value in a sorted vector,16Sign Value> encoded as 2≥ encoded as 1= encoded as 0≤ encoded as -1< encoded as -2Table 4.1.2: Sign value encoding for inequalities for bounded and unbounded regionsreturns the index of location where the value is) was found [9]. As this algorithm onlysupports vectors, the method was expanded to support a row search in a matrix, such thatthe algorithm returns the index of a row match for a given inequality. Utilizing this method,the two models can be merged and for duplicate entries, as the binary search returns theindex of the initial duplicate row, the reference is updated such that the model being mergedhas its duplicated entries updated with the indexes from the base model.In looking forward to other parts of the algorithm, additional information needs to beencoded during the conversion algorithm. Unlike the Voronoi diagram, the Delaunay tes-sellation does not maintain uniqueness in terms of using a center to treat the problem as apoint location problem. In order to accelerate the evaluation, additional methods have beenintroduced such that each triangulation also has an associated center of gravity or centroidas seen in Model (2.2.3). Thus, the centroid is computed for each face and stored in thetriangulation row entry.Additionally, a neighbour list is constructed during the conversion that allows each tri-angulation to maintain a list of other triangulations it is next to. For each triangulation ti,an index list of neighbouring triangulations is maintained such that the neighbour will onlyappear in the list if and only if it has a common edge with ti. The centroid and neighbourlist methods have not been fully implemented for the unbounded regions.17Algorithm 5: PLQ2 convert to IPLQ2 AlgorithmInput: PLQ2Output: as inequality-based PLQ2begininitialization;forall triangulations ti ∈ D in PLQ2 doCompute the linear-equalities that form the triangulation;Compute the centroid for the face;forall edges in the triangulation dorecord the index of the neighbouring triangulation;endCheck for duplicate entries in the existing list of inequalities;endCompute the convex hull on P forming k + 1 hull points and k edges;forall edges k in the convex hull docompute the 2 bounding inequalities and build 2k unbounded faces with atmost 2k inequalities;endforall unbounded inequalities docheck to see if it exists in the list of bounded equalities and if it does updatethe index in the face inequality reference list;endConcatenate the two models;Return the inequality-based IPLQ2;endProposition 4.1.1. The PLQ2 convert to IPLQ2 algorithm runs in O(n2).We now recall some factsa136 Let P be a set of n points in the plane, not all collinear, and let k denote the numberof points in P that lie on the boundary of the convex hull of P. Then any triangulationof P has 2n−2−k triangles and 3n−3−k edges.a136 The Delaunay triangulation contains at most O(nd/2) simplexes, thus ina822, for n inputpoints, there will be at most O(n) triangulations. For d = 2 the model will have atmost O(k) triangulations.a136 In the plane a822, if there are b vertices on the convex hull, then any triangulation ofthe points has at most 2n−2−b triangles, plus one exterior face [8].a136 The centroids of a set of n points in d dimensions can be computed trivially in O(dn) =O(n) [12, p. 106].We are now ready to prove Proposition 4.1.118Proof. Let Q be a bounded bivariate PLQ model with n input points with at most n trian-gulations ti ∈ D. For each ti, there are three vertices that form the face thus the face hasthree edges for which three linear-inequalities are needed to represent the triangulation. Theequation of each edge is computed in linear time using the point slope and intercept form foreach set of points. Each face can be represented by 3 inequalities in O(1) and all inequalitiescalculated are n∗O(1) = O(n). Checking to see if a single inequality exist in the inequalitylist is in O(logm) where m is the number of inequalities. It follows, as the model has at mostn triangulation, and each has 3 edges, that m = 3n and O(logm) = O(log(3n)) = O(logn);thus, for n faces, checking for duplicates is O(nlogn).To compute the centroid of the new triangulation, the centroid for a single triangulationcan be computed in O(p) where p = 3 and is the number of vertices in the triangulation.Thus, for n triangulations, the complexity is O(3n) = O(n).To compute the neighbour list, each triangulation maintains three vertices. For each setof vertex that participates in an edge, we index into the triangulation matrix and recoverthe members in at worst O(n). A bottleneck is introduced by the fact that a single vertexcould be a member of every triangulation in the model; hence its triangulation membershipis n. This cycle is repeated for each vertex in the triangulation; thus for n triangulations,the worst case complexity for building the neighbour list is O(n∗n) = O(n2).To compute the unbounded regions, by Definition 1.1.8 and 1.1.5 the Delaunay tessella-tion is a maximal planar subdivision, thus each point in a822 is a vertex. It then follows thatthere are b vertices and b number of points, thus b = n and the number of triangles is2n−2−b + 1 = 2n−2−n + 1= n−2 + 1= n−1.Thus, the number of triangles is at most n−1 for D. Let k be the number of edges in theconvex hull of D. The number of edges in the convex hull is be related to the number oftriangles as 2m − 2 − k where m is the number of input points. Thus, for the unboundedPLQ where n is the number of triangles, the number of edges is n = 2m−2−k. With m = nwe have n = 2n−2−k so k = n−2 = O(n).There are at most n−2 edges in the convex hull of D and for each edge, two inequalitiesare required which can be each computed in O(1). Hence, at most 2k = 2n−4 inequalitiesare required, each taking O(1), and the inequalities for the unbounded region are computedin O(n) with at most 2n−4 inequalities.In joining the two models, the existence of an unbounded inequality in bounded equalitylist can be checked in O(logn). Thus 2n−4 inequalities can be checked in O((2n−4)logn) =O(nlogn). With any duplicate inequalities removed, the two models can be concatenated inO(1). It follows that the run time complexity for the conversion of a bounded tessellation-based bivariate PLQ to an unbounded inequality-based bivariate PLQ without the neighbourlist structure or centroid list is O(n) + O(nlogn) + O(n) + O(nlogn) + O(1) = O(nlogn).Additionally, the run time complexity for the conversion of a bounded tessellation-basedbivariate PLQ to an unbounded inequality-based bivariate PLQ with the neighbour list19structure or centroid list is O(n)+O(n)+O(n2)+O(nlogn)+O(n)+O(nlogn)+O(1) =O(n2).4.2 Evaluation AlgorithmThe evaluation of an unbounded inequality-based bivariate PLQ over a set of points (xi,yi) ∈X is achieved by searching through the triangulations formed by the inequalities and deter-mining if a specific point is within a given region. MATLAB does not support functionalityto search through triangulations as defined by inequalities, thus for the initial method, thetriangulations will be toured through in order. Fortunately, the triangulations have an in-herent natural ordering from the build function. By definition, if a point (x,y) falls on theboundary of a triangulation, it can be evaluated as part of either facet. If a point thatis being evaluated is not in the domain of the model, the index is returned as +∞. Twodifferent methods are presented. Only Method A will be proven as Method B is not fullyimplemented.In order to accelerate the evaluation, points are ordered in lexicographic order to exploitthe fact that in the process of the build operation, there is an inherent ordering of trian-gulations. Thus, it is possible, if both X and T share the same ordering, to determine thetriangulation in which each point is located in a linear fashion. Once a triangulation is foundto bound the point under evaluation, it is evaluated using the associated linear-quadraticfunction.Algorithm 6: IPLQ2 eval AlgorithmInput: IPLQ2,XOutput: a vector of values Zbegininitialization;Order X in increasing lexicographic order;Set i,j = 0;for each xi ∈ X doif xi satisfies the constraints of the triangulation tj thenevaluate with associated linear-quadratic to Z and check xi+1;endcheck xi in tj+1;if j ≥ size(T ), the number of triangulations thenreset j = 0;endendReorder Z to original order;Return Z;endA second approach was investigated to improve the performance. Using the centroidsfor each face, it was conjectured that the problem could be framed as a point location20type problem. Even though locating a point to a given center may not absolutely locate apoint within a face, we conjecture that this significantly decrease the search space. Oncea point has been located to a centroid, the point is evaluated under the bounds of thetriangulation associated with the given centroid. If the point is not in the bounded region,we take the next closest centroid and so on. The points under evaluation are still ordered inincreasing lexicographic order, and once a triangulation is found that contains a point, thenall subsequent points are evaluated under that triangulation until a point is found that isnot in the bounded face. At this point, the centroid search is computed again and repeateduntil all points have been evaluated.Algorithm 7: IPLQ2 eval v2 AlgorithmData: IPLQ2,XResult: a vector of values Zbegininitialization;Order X in increasing lexicographic order;Set i,j = 0;forall xi ∈ X doSearch the j triangulations for the closest centroid;if the point is in triangulation t thenevaluate with associated linear-quadratic to Z;check xi+1 is in the face;endcheck xi in the next closest face;endReorder Z to original order;Return Z;endWe first recall known resultsa136 Sorting a set of items has a computational run-time complexity of O(nlogn) and astorage complexity of O(n).a136 Searching a set of sorted items has a computational run-time complexity of O(logn).a136 If a point p is inside the convex hull of a Delaunay triangulation, its nearest vertexneed not be one of the vertices of the triangle containing p.Proposition 4.2.1. The evaluation of k points in a model of n triangulations using Algo-rithm 6 takes O(kn).Proof. Let k be the number of points (xi,yi) ∈ X where 0 ≤ i < k to be evaluated ina bivariate inequality-based model IPLQ2, that contains n triangulations. Ordering of Xis accomplished using a standard sorting algorithm which runs in O(klogk). For a single(xi,yi), the worst case is that n triangulations must be searched until a legitimate bounded21triangulation can be found. If this is the case for every xi ∈ X, then locating the validtriangulation for each point is O(kn).Once a legitimate bounded triangulation is found for a point, the linear-quadratic functioncan be evaluated in O(1) as stated in Proposition 3.2.1. It then follows that the evaluationof k points can be completed in O(k) which is also a tight bound as all k points mustbe evaluated; the linear-quadratic function evaluation is Θ(k). Returning the points totheir original order is O(k). The worst-case run time complexity for evaluating k points isO(klogk) + O(kn) + O(8) + O(k) = O(kn).Conjecture 4.2.1. Algorithm 7 has a run-time complexity of O(n)4.3 Addition4.3.1 AlgorithmThe IPLQ2 add algorithm takes as input, two unbounded inequality-based bivariate PLQmodels of size n and returns an inequality-based bivariate PLQ that is the addition of thetwo input models. The addition method for the inequality-based model is functionally verysimilar to the bounded addition but operationally it is quite different. The addition algorithmstarts with one inequality-based PLQ as a base model PLQ2a. To PLQ2a, the inequalitiesfrom the second PLQ, PLQ2b are added one-by-one. As each inequality is added, it willintersect n triangulations and divide into 2n convex polytope tessellations. The resultingsplit may not be triangular and on each new tessellation, the centroid and neighbour list iscalculated. Once all of the inequalities have been added, the resulting model contains theunion of the tessellations from PLQ2a and PLQ2b, but the linear-quadratic function valuefor each tessellation must be computed.Using the new centroid list Cprime, the centroids are evaluated under PLQ2a and PLQ2bexcept instead of evaluating the point, the linear-quadratic function is extracted from eachand the summation of the two linear-quadratic equations stored in the relative tessellationassociated with each ci ∈ Cprime. We conjecture that this problem can be formulated as a mapoverlay problem which Becker et al. have shown to have a complexity of Θ(n2)[5]. Thismethod is not fully implemented at this time.The addition method utilizes the neighbour list to determine the direction of intersectionwhen an inequality is added to streamline the addition process. The algorithm is as follows:22Algorithm 8: IPLQ2 add AlgorithmInput: IPLQ2a,IPLQ2bOutput: IPLQ2 which is the summation of the two input IPLQ2begininitialization;Order X in increasing lexicographic order;Set i,j = 0;forall inequalities ii ∈ I in PLQ2b doAdd the inequality to PLQ2a as PLQ2c;Locate a triangulation cut by the inequality;forall triangulations bisected by the inequality doform two new tessellations for each face split;Compute the centroids for each new face formed;endendEvaluate the centroids for each face under PLQ2a and extract the linear-quadraticfunctions;Evaluate the centroids for each face under PLQ2b and extract the linear-quadraticfunctions;Sum the two sets of linear-quadratic functions;Update the linear-quadratic functions for each tessellation in PLQ2c;Return PLQ2c;endProposition 4.3.1. The addition algorithm has a computational run time complexity ofO(n3) and a computational storage complexity of O(n2).Proof. Let PLQ2a be an unbounded inequality-based bivariate PLQ model of n triangula-tions and let PLQ2b be an unbounded inequality-based bivariate PLQ model of m trian-gulations. As PLQb have n triangulations, by conjecture 4.1.1, it contains no more than2n−4 = O(n) inequalities. For each inequality in PLQb, it is added to PLQa which containsn triangulations and if n triangulations are bisected, 2n tessellations are formed from the23addition of a single inequality. Thus adding n inequalities forms at most n2 new facesn +2nsummationdisplayi=n+1i = n +nsummationdisplayi=1(i + n)= n +nsummationdisplayi=1i +nsummationdisplayi=1n= n + (n + 1)(n)2 + n2= n + n22 +n2 + n2= 3n22 +3n2 = O(n2)For each inequality that is added, 2n tessellation are formed and 2n centroids are com-puted in O(n) as stated in Section 5; thus for n inequalities added, n2 centroids Cprime arecomputed. To extract the linear-quadratic function from PLQ2a and PLQ2b, C is eval-uated under each model, returning the function instead of the resultant of the function.By Proposition 4.3.1 the extraction of the linear-quadratics functions is in O(kn), wherek = n2 thus the complexity is O(2n3). Updating the linear-quadratic function in PLQ2ccan be completed in O(n2) as it is the operation of indexing into a matrix of size n2.Thus, the complete complexity for adding two unbounded inequality-based bivariate PLQ’sis O(n2) + O(n3) + O(n2) = O(n3).Corollary 4.3.1. If the complexity of IPLQ2 eval is O(n+k) then the complexity of IPLQ2 addis O(n2).4.4 Scalar Multiplication AlgorithmScalar multiplication for the inequality-based PLQ is a straightforward linear-time algorithmin which the linear-quadratic equation for each face in the PLQ model in multiplied by ascalar value. It does not effect any other component of the data structure.24Algorithm 9: IPLQ2 mult AlgorithmInput: IPLQ2,cOutput: an IPLQ2 that has been scalar multiplied by cbegininitialization;forall each triangulation ti ∈ D in PLQ2 domultiply the linear-quadratic equation at index i by the constant c;endReturn PLQ2;endProposition 4.4.1. The IPLQ2 mult algorithm runs in O(n).25Chapter 5Conclusions and Future Work5.1 ConclusionsThe following complexities were found for each model:OperationModel build eval add multTessellation O(nlogn) O(nk) O(n2) O(n)Inequality O(n2) O(nk) O(n3) O(n)Voronoi1 O(nlogn) O(n + k)2 O(n2) O(n)Table 5.1.1: Complexity Results for Fundamental Convex TransformsIt was observed empirically that, for functions that contained large amounts of linear datathe algorithms tended to perform better than the upper bound. Overall, for linear encodeddata, the models perform empirically better than the worst case complexity. Encodingstrictlyconvexfunctionsdoespresentachallengeforthis modelas little ornosimplificationorreduction in the size of the input compared to the output can occur, leading to a performancebottleneck.It was also observed that most methods that generate or permeate the model are outputsensitive in terms of the number of faces generated. As previously mentioned, as this is azero order model, strictly convex data cannot be efficiently encoded. Thus, as the modelimproves and is able to encode quadratic information, a performance increase may be seenas the number of faces required to encode data over a given region may decrease.Limitations also exist with the tessellation-based model for strictly convex functionswith a large number of parts in a small area. The convex hull algorithm starts to reach aperformance barrier when the angles of each triangulation become extremely small, such thattrouble is encountered when trying to resolve the location of a point. With the introductionof a model that allows for the encoding of quadratic information, this should no longer be abarrier to performance. A possible way to extend the tessellation-based model to functions1Results for the Voronoi model are conjectured.2Since both the Vertices and the list of points to evaluate can are sorted, this amounts to merging tosorted sequences (n+k) + sorting : n log n + k log k but the vertices can be maintained sorted i.e. sortedin the build. So the cost drops to n+ k log k. If the set of k points is a grid (as is most often the case), it isassumed already sorted along each axis (so trivially sorted in the plane) so the cost is now n+k.If nothing is sorted, evaluation can still be done in k log n by repeatedly searching in the sorted list ofvertices.26with unbounded domains would be to store vertices and half-lines separating the unboundedfaces (e.g. each vertex is associated with a list of directions corresponding to the edges ofthe tessellation going through that vertex). Such information is redundant for vertices inthe interior of the domain (since it can be extracted from the other vertices) but is criticalfor vertices on the boundary of the domain to handle unbounded faces.AdditionalworkneedstobecarriedouttocompletetheimplementationfortheIPLQ2 evaland IPLQ2 add methods as well as verifying and improving the run time complexity. Weknow that the addition algorithm when formulated as a map overlay type problem can havea run time complexity of O(n2), while our method is still bounded by O(n3) with the bot-tleneck being the extraction and summation of the linear-quadratic equation for each newface. Future Work will focus on streamlining and improving this functionality to bring themethod within the Ω(n2) bound. Further work will also be conducted into improving theevaluation method with improved search methods.Further work will investigate if a PLQ function can be encoded using a Voronoi diagram,which will be computed from the initial Delaunay tessellation. Voronoi diagrams have previ-ously been used to encode a convex hull as a finite set of intersecting half-spaces [3]. A userwould define the point value in the model and the algorithm would compute and project aVoronoi diagram into a822. Encoded with each center would be the value of the function asa linear-quadratic equation. Thus a PLQ function in a823 is encoded as a matrix containingVoronoi centers (xi,yi) ∈ V in a822 associated with the corresponding linear-quadratic equa-tion. Equation (2.1.1) encodes the function value of each part of the PLQ to the boundariesof the polytope for the Voronoi diagram in a822 containing the Voronoi centers V.Adding the location of each Voronoi center (xi,yi) ∈ V with Equation (2.1.2) a PLQfunction is encoded into a single matrix asx1 y1 a1,1 a1,2 a1,3 a1,4 b1,1 b1,2 c1x2 y2 a2,1 a2,2 a2,3 a2,4 b2,1 b2,2 c2· · · · · · · ·· · · · · · · ·xi yi ai,1 ai,2 ai,3 ai,4 bi,5 bi,6 ci. (5.1.1)This model has desirable characteristics due to the properties of Voronoi diagrams. Aspreviously stated, the Voronoi diagram of a set on n sites in the plane can be computed witha line sweep algorithm in O(nlogn) and with a storage complexity of O(n) [8, 3] with thestorage complexity for the complete model being O(n).The desirability of this model is further enhanced in terms of evaluating the bivariatePLQ over a set of points, xi,yi ∈ X. This problem can be treated as a variant of the post-office problem [8, 3] with solutions for a single point query in O(nlogn) [3], noting that witha Voronoi diagram, each Voronoi cell is uniquely defined by xi,yi ∈ V, and by definition, agiven cell contains all points closer to va than to any other center vi ∈ V.We feel that this model may provide a good basis for a bivariate PLQ in a823, but uncer-tainty exists regarding the addition function. We anticipate that the addition of two Voronoimodels may be treated as a map overlay problem and future work will focus on investigation27of this function for this possible model. Provided that this is the case, then we will have twopossible models to represent bivariate PLQ functions.As the model presented only represents zero order data, the data structure can becomelarge for strictly convex functions. With the introduction of a first and second order model,improved encoding will occur to reduce the storage complexity. Additional limitations existin the current model for dense tessellation where the linear-quadratic function value betweentriangulations is very small, such that they can be joined into a single face to reduce theoverall complexity of the model; thus, a clean function needs to be introduced to allow forthe accommodation of such functionality.Questions still exist in terms of the ratio between the input size k and the output size nand further work needs to be done to investigate this relationship. With respect to the codelibrary, continual improvement is needed to add support for degenerate cases in additionto improving test coverage and unit testing. Additionally, support for other fundamentalconvex transforms needs to be expanded. Work to parallelize methods may be undertakento improve the run time performance of addition and conversion methods.28Chapter 6APPENDIXNumerical Examples The following sections include numerous examples of using the specificmethods previously introduced. All examples are run in MATLAB R2007b.6.1 Building a Tessellation-Based PLQ fromPointwise Data6.1.1 ExampleBuilding a tessellation-based PLQ from pointwise data involves two steps. First, the functiondata must be formed as two separate elements; the first, X, is a 2xn matrix that representsthe points in a822 that the function is defined over. The point sampling does not have tobe regular. The second element is a vector f(X) which represents the function value overX. Building the function involves a call to PLQ2 build with X and f(X) being passes asparameters. The function will return a structure the contains two elements. The first elementis a matrix X, the vertices for the triangulations that define the PLQ. The second elementis a matrix TES that contains the vertex indices for each triangulation and the linear-quadratic function associated with each triangulation. It should be noted that MATLABuses 1 based indexing for arrays and matrices. A helper method PLQ2 show tes whichtakes as an argument a tessellation-based PLQ2, will produce a graph is a822 that graphicallyrepresents that triangulation produced.%define a series of points and the function value>> [x y] = meshgrid(-10:2:10);>> x = reshape(x,size(x,1)^2,1);>> y = reshape(y,size(y,1)^2,1);>> f = abs(x) + abs(y);>> X = [x y];>> plq2 = PLQ2_build(X ,f);%The structure of the model>> plq2plq2 =TES: [8x10 double]X: [9x2 double]%The vertexes>> plq2.X29(a) Tessellation ina822 produced by the build method (b) Run Times for the PLQ2 build MethodFigure 6.1.1: Build Method Resultsans =-10 -10-10 0-10 100 -100 00 1010 -1010 010 10%The faces>> plq2.TESans =2 4 5 0 0 0 0 -1 -1 01 4 2 0 0 0 0 -1 -1 05 4 8 0 0 0 0 1 -1 08 4 7 0 0 0 0 1 -1 06 2 5 0 0 0 0 -1 1 03 2 6 0 0 0 0 -1 1 05 8 6 0 0 0 0 1 1 06 8 9 0 0 0 0 1 1 0>> PLQ2_show_tes(plq2);The results of the call to the helper function PLQ2 show tes can be seen in figure 6.1(a).30Figure 6.1.2: Evaluation of a tessellation-based PLQ26.1.2 Commentsinput size 25 441 1681 6561 40401 160801number of triangulations 32 800 3200 12800 80000 320000time(seconds) 0.0071 0.1943 0.4179 1.6626 9.8149 61.9158Table 6.1.1: Complexity Results of the Build Method for a Strictly Convex FunctionIn looking at the empirical run time complexity from timing results from MATLAB as inTable 6.1.1 and in Figure 6.1(b), the function appears to progress beneath the conjecturedrun time complexity of O(nlogn). The storage complexity also appears to be in O(n) as thenumber of triangulations generated increases linearly with the number of input points.6.2 Evaluating a Bounded PLQ6.2.1 ExampleEvaluating tessellation-based PLQ from pointwise data involves two steps. First, a matrixof the points X in 2 x n over which the function is to be evaluated. Once the points X havebeen defined, the second step is to invoke the function involves with call to PLQ2 eval withPLQ2 and X being passed as parameters. The parameters PLQ2 is the model over whichX is being evaluated. The function will return a vector of function values for each xi ∈ X%using the model PLQ2 defined previously%define a matrix of points to be evaluated>> ti = -10:10:10;>> [XI,YI] = meshgrid(ti,ti);>> xx = reshape(XI,size(XI,1)^2,1);>> yy = reshape(YI,size(YI,1)^2,1);31>> XX = [xx yy];%and evaluate>> z = PLQ2_eval(plq2,XX);>> zz =20102010010201020To plot the results graphically, a regular grid can be defined as above and plotted as>> ZI = griddata(xx,yy,z,XI,YI);>> mesh(XI,YI,ZI), hold>> g2 = plot3(xx,yy,z,’.’),hold offCurrent plot heldg2 =160.0221with the results of the plot show in figure 6.1.1.6.2.2 CommentsThe empirical results have previously been presented in Section 3.2 and it was seen that theevaluation is bounded by O(nk)6.3 Adding Two Bounded PLQ’s6.3.1 ExampleAdding two tessellation-based PLQ’s is accomplished with call to PLQ2 add with two PLQ’s,PLQa and PLQb being passed as parameters. The function will return a PLQ that is thesummation of the two models. The two models do not have to be defined over the samedomain.%define model>> [x y] = meshgrid(-10:5:10);>> x = reshape(x,size(x,1)^2,1);>> y = reshape(y,size(y,1)^2,1);>> f = abs(x) + abs(y);>> X = [x y];32>> plq2a = PLQ2_build(X ,f);>> size(plq2a.TES)ans =8 10%define model 2>> [x y] = meshgrid(-10:5:10);>> x = reshape(x,size(x,1)^2,1);>> y = reshape(y,size(y,1)^2,1);>> f = x.^2 + y.^2;>> X = [x y];>> plq2b = PLQ2_build(X ,f);>> size(plq2b.TES)ans =32 10%and add the two models>>plq2c = PLQ2_add(plq2a,plq2b);>>size(plq2c.TES)ans =32 106.3.2 CommentsIn this example, two functions are added producing a third. Due to the increasing size of thedata structure, the results are show graphically. Figure 6.1(a) and Figure 6.1(b) representthe two functions that are being added together. The results of the summation are shownin Figure 6.1(c).An important observation must be made with the bounded tessellation-based PLQ2 add;if the two models have different domains, the resulting domain will be the intersection of thetwo models as the model only represents values that are defined. Hence, it is possible to havea reduction in number of faces in the model due to this property. Another observation wasthat in testing, the bounded tessellation-based PLQ2 add tended to run in O(n) as in Figure6.1(d) suggesting that the average run time complexity may be less that the theoretical worstcase run time complexity of O(n2).33(a) PLQa (b) PLQb(c) PLQc: the Results of the Addition of PLQa andPLQb(d) Run Time for the Addition of Two PLQ ModelsFigure 6.3.1: The Addition of Two PLQ Models346.4 Scalar Multiplication of a tessellation-based PLQ6.4.1 ExampleMultiplyingatessellation-basedPLQbyascalarconstantiscomputedbyacalltoPLQ2 multwith PLQ2 and k being passes as parameters where k is the scalar constant. The functionwill return a tessellation-based PLQ.>> plq2plq2 =TES: [8x10 double]X: [9x2 double]>> plq2.TESans =2 4 5 0 0 0 0 -1 -1 01 4 2 0 0 0 0 -1 -1 05 4 8 0 0 0 0 1 -1 08 4 7 0 0 0 0 1 -1 06 2 5 0 0 0 0 -1 1 03 2 6 0 0 0 0 -1 1 05 8 6 0 0 0 0 1 1 06 8 9 0 0 0 0 1 1 0>> plq2 = PLQ2_mult(plq2,10);>> plq2.TESans =2 4 5 0 0 0 0 -10 -10 01 4 2 0 0 0 0 -10 -10 05 4 8 0 0 0 0 10 -10 08 4 7 0 0 0 0 10 -10 06 2 5 0 0 0 0 -10 10 03 2 6 0 0 0 0 -10 10 05 8 6 0 0 0 0 10 10 06 8 9 0 0 0 0 10 10 06.4.2 CommentsVisually, the results can seen in figure 6.1(a) and 6.1(b) where a function has been multipliedby a constant, k = 10. In looking at the empirical run time complexity as in figure 6.1(c) isseen to be linear with an increasing number of triangulations in the model. This supportsconjecture 3.4.1 that the algorithm runs in O(n) in terms of the number of triangulation inthe model.35(a) PLQ Before Multiplication by a Scalar Constant (b) PLQ After Multiplication by a Scalar Constant(c) PLQ2 Mult Run Time ComplexityFigure 6.4.1: PLQ Multiplication36Figure 6.5.1: Run Time Complexity for Model Conversion6.5 Converting a tessellation-based PLQ to aninequality-based PLQ6.5.1 ExampleConverting from a bounded tessellation-based PLQ model to an inequality-based PLQ modelcan be done with a single call to the method PLQ2 convert to IPLQ2 which will return aninequality-based PLQ. In the example, a simple tessellation-based PLQ is converted to aninequality-based PLQ.>>plq2 =TES: [8x10 double]X: [9x2 double]>> iplq2 = PLQ2_convert_to_IPLQ2(plq2);>> iplq2iplq2 =ineq: [20x3 double]faces: {12x2 cell}It should be observed that the number of faces in the model has increased from 8 to 12 forthis example. This is due to the fact that the tessellation-based model is a bounded modeland when converted to an inequality-based model, the domain is expanded such that theentire a822 space is encoded; hence the need for additional faces.376.5.2 CommentsA simulation was run over a series of models with an increasing number of tessellations upto 20,000 faces. The model that was used was a strictly convex function such that it wouldgenerate a worst case scenario. The results are seen in figure 6.5 which shows processortime plotted against the size of the model. Additionally, the theoretical O(n2) is shown asa comparison. It appears that the empirical run time complexity supports Proposition 4.1.1such that the worst case run time complexity is O(n2).386.6 Evaluating an inequality-based PLQ6.6.1 ExampleEvaluating tessellation-based IPLQ from pointwise data involves two steps. First, a matrixof the points X in 2 x n over which the function is to be evaluated. Once the points X havebeen defined, the second step is to invoke the function involves with call to IPLQ2 eval withIPLQ2 and X being passes as parameters. The parameter ILQ2 is the model over whichX is being evaluated. The function will return a vector of function values for each xi ∈ X%using the model IPLQ2 defined previously>> ti = -10:10:10;>> [XI,YI] = meshgrid(ti,ti);>> xx = reshape(XI,size(XI,1)^2,1);>> yy = reshape(YI,size(YI,1)^2,1);>> XX = [xx yy];%and evaluate>> z = IPLQ2_eval(iplq2,XX);>> zz =201020100102010206.6.2 CommentsTwo separate tests were run, one looking at the empirical run time complexity with a varyingnumber of input points over a fixed number of tessellation. The second evaluated the runtime complexity with a fixed number of points over a varying number of tessellations. Bothexhibited linear behavior over the scope of the test. A third test was run looking at varyingboth the number of points k and the number of tessellations k. Empirically, the inequal-ity evaluation performs in a linear fashion over nk as seen in figure 6.1(c) which supportsconjecture 4.2.16.7 Additional MethodsUsage of the additional methods supporting multiplication and addition for the inequalitymodel follow the examples for the tessellation-based model.39(a) Evaluation with a Constant Number of Points (b) Evaluation with a Constant Number of Faces(c) IPLQ2 eval Run Time Complexity Over nkFigure 6.6.1: IPLQ Evaluation40Bibliography[1] Matlab r2007b, 2008. http://www.mathworks.com/.[2] Scilab, 2008. http://www.scilab.org/.[3] F. Aurenhammer, Voronoi diagrams - a survey of a fundamental geometric datastructure, ACM Comput. Surv., 23 (1991), pp. 345–405.[4] C. B. Barber, D. P. Dobkin, and H. Huhdanpaa, The quickhull algorithm forconvex hulls, ACM Trans. Math. Softw., 22 (1996), pp. 469–483.[5] L. Becker, A. Giesen, K. Hinrichs, and J. Vahrenhold, Algorithms for per-forming polygonal map overlay and spatial join on massive data sets, ”Lecture Notes inComputer Science”, ”1651” (1999), pp. 270–285.[6] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press,Cambridge; United Kindgdom, 2006, c2004.[7] T. M. Chan, Optimal output-sensitive convex hull algorithms in two and three dimen-sions, Discrete Comput. Geom., 16 (1996), pp. 361–368. Eleventh Annual Symposiumon Computational Geometry (Vancouver, BC, 1995).[8] M. de Berg, M. van Kreveld, M. Overmars, and O. Schwarzkopf, Com-putational Geometry: Algorithms and Applications, Springer-Verlag Berlin Heidelberg,New York; New York, second ed., c2000, c1997.[9] M. Khan, Binary search for numeric vector, 2006.http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=11287.[10] A. J. King and D. L. Jensen, Linear-quadratic efficient frontiers for portfoliooptimization, Applied Stochastic Models and Data Analysis, 8 (1992), pp. 195–207.http://dx.doi.org/10.1002/asm.3150080309.[11] Y. Lucet, H. H. Bauschke, and M. Trienis, The piecewise linear-quadratic modelfor computational convex analysis, Comput. Optim. Appl., (2007). Online publication.[12] F. P. Preparata and M. I. Shamos, Computational Geometry, Springer-VerlagBerlin Heidelberg, New York; New York, c1985.[13] A. Rantzer and M. Johansson, Piecewise linear quadratic optimal control, 1997.citeseer.ist.psu.edu/rantzer97piecewise.html.41
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Undergraduate Research /
- The Piecewise Linear-Quadratic Model for Convex Bivariate...
Open Collections
UBC Undergraduate Research
The Piecewise Linear-Quadratic Model for Convex Bivariate Functions Fazackerley, Scott Ronald 2008
pdf
Page Metadata
Item Metadata
Title | The Piecewise Linear-Quadratic Model for Convex Bivariate Functions |
Creator |
Fazackerley, Scott Ronald |
Date Issued | 2008 |
Description | A Piecewise Linear-Quadratic (PLQ) function is used to describe data that can be represented by continuous functions with a piecewise linear domain for which the function is either linear or quadratic on, each piece of its domain [11]. While extensive work has been done with PLQ models for convex univariate functions, my work investigates the development of a two dimensional model that allows the implementation of algorithms for computing fundamental convex transforms. PLQ functions have been described in the literature, the efficient implementation of the algorithms requires the careful selection of a data structure. Initial investigation examined using a Voronoi diagram to represent the projection of the bivariate function into R2, to take advantage of efficient algorithms for the point location problem [3]. While suitable for representing a single PLQ, with operations for building a zero order model from data and evaluating the function over an irregular grid, we are uncertain if it can be extended to compute other fundamental convex transforms efficiently. In examining the Voronoi model, we were able to extend the base concept to represent the bivariate PLQ using two different representations: a tessellation-based model and a linear inequality-based model. The tessellation model is restricted to representing a PLQ with a bounded domain. The model represents data through triangular faces in a tessellation in R2 where each face is defined by its vertices and the associated function value. This model allows for the efficient evaluation over an irregular grid, the addition of two PLQ functions with bounded domains and multiplication by a scalar, but does not allow for the representation of an unbounded domain. To allow for an unbounded domain, a dual model was developed using linear inequalities to represent each face in R2. Numerical results are presented for each model and computational complexity of model components are discussed. |
Genre |
Other |
Type |
Text |
Language | eng |
Series | University of British Columbia. COSC 449 |
Date Available | 2010-07-30 |
Provider | Vancouver : University of British Columbia Library |
DOI | 10.14288/1.0052224 |
URI | http://hdl.handle.net/2429/27051 |
Affiliation |
Irving K. Barber School of Arts and Sciences (Okanagan) |
Peer Review Status | Unreviewed |
Scholarly Level | Undergraduate |
Aggregated Source Repository | DSpace |
Download
- Media
- [if-you-see-this-DO-NOT-CLICK]
- Thesis Scott Fazackerley.pdf [ 554.46kB ]
- Metadata
- JSON: 1.0052224.json
- JSON-LD: 1.0052224+ld.json
- RDF/XML (Pretty): 1.0052224.xml
- RDF/JSON: 1.0052224+rdf.json
- Turtle: 1.0052224+rdf-turtle.txt
- N-Triples: 1.0052224+rdf-ntriples.txt
- Original Record: 1.0052224 +original-record.json
- Full Text
- 1.0052224.txt
- Citation
- 1.0052224.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
China | 3 | 0 |
Canada | 2 | 0 |
Germany | 1 | 0 |
United States | 1 | 0 |
City | Views | Downloads |
---|---|---|
Shenzhen | 3 | 0 |
Montreal | 2 | 0 |
Unknown | 1 | 0 |
Ashburn | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.52966.1-0052224/manifest