Arts and Sciences, Irving K. Barber School of (Okanagan)
Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)
DSpace
UBCO
Kumar, Deepak
Master of Science - MSc
University of British Columbia
Computing the convex envelope or biconjugate is the core operation that bridges the domain of nonconvex analysis with convex analysis. For a bivariate PLQ function de ned over a polytope, we start with computing the convex envelope of each piece. This convex envelope is characterized by a polyhedral subdivision such that over each member of the subdivision, it has an implicitly de ned rational form(square of a linear function over a linear function). Computing this convex envelope involves solving an exponential number of subproblems which, in turn, leads to an exponential time algorithm. After that, we compute the conjugate of each such rational function de- ned over a polytope. It is observed that the conjugate has a parabolic subdivision such that over each member of its subdivision, it has an implicitly de ned fractional form(linear function over the square root of a linear function). This computation of the conjugate is performed with a worst-case linear time complexity algorithm. Finally, some directions and insights about computing the maximum of all the conjugates of each piece of a PLQ function and then the conjugate of that to obtain the biconjugate are provided as conjectures for future work.
Conjugate of some Rational functionsand Convex Envelope of Quadraticfunctions over a PolytopebyDeepak KumarBachelor of Technology, Indian Institute of Information Technology, Design andManufacturing, Jabalpur, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Computer Science)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)January 2019c© Deepak Kumar, 2019The following individuals certify that they have read, and recommendto the College of Graduate Studies for acceptance, a thesis/dissertation en-titled:Conjugate of some Rational functions and Convex Envelopeof Quadratic functions over a Polytopesubmitted by Deepak Kumar in partial fulfilment of the requirements ofthe degree of Master of Science.Dr. Yves Lucet, I. K. Barber School of Arts & SciencesSupervisorDr. Heinz Bauschke, I. K. Barber School of Arts & SciencesSupervisory Committee MemberDr. Warren Hare, I. K. Barber School of Arts & SciencesSupervisory Committee MemberDr. Homayoun Najjaran, School of EngineeringUniversity ExamineriiAbstractComputing the convex envelope or biconjugate is the core operation thatbridges the domain of nonconvex analysis with convex analysis. For a bi-variate PLQ function defined over a polytope, we start with computing theconvex envelope of each piece. This convex envelope is characterized by apolyhedral subdivision such that over each member of the subdivision, ithas an implicitly defined rational form(square of a linear function over alinear function). Computing this convex envelope involves solving an expo-nential number of subproblems which, in turn, leads to an exponential timealgorithm.After that, we compute the conjugate of each such rational function de-fined over a polytope. It is observed that the conjugate has a parabolicsubdivision such that over each member of its subdivision, it has an im-plicitly defined fractional form(linear function over square root of a linearfunction). This computation of the conjugate is performed with a worst-caselinear time complexity algorithm.Finally, some directions and insights about computing the maximum ofall the conjugates of each piece of a PLQ function and then the conjugate ofthat to obtain the biconjugate are provided as conjectures for future work.iiiLay SummaryOptimization problems occur in a wide variety of fields ranging fromEngineering to Mathematical Finance and solving these problems play acrucial role in the process. When the optimization problem is nonconvex itis often difficult to solve. In this thesis, we present a method to convert someclass of nonconvex problems to convex problems, which are often easier tosolve, and a method to compute their conjugates which can help solving theproblem faster than before.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Notations and Abbreviations . . . . . . . . . . . . . . . xiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 PLQ Functions . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Convex Envelope . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Fenchel conjugate . . . . . . . . . . . . . . . . . . . . . . . . 41.4 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5Chapter 2: Preliminaries and Notations . . . . . . . . . . . . . 8Chapter 3: Convex Envelope of Bivariate Quadratic Func-tions over a polytope . . . . . . . . . . . . . . . . . 133.1 Quadratic Reduction . . . . . . . . . . . . . . . . . . . . . . . 133.2 Problem Formulation[Loc16] . . . . . . . . . . . . . . . . . . . 163.2.1 Problem Formulation . . . . . . . . . . . . . . . . . . 163.2.2 Expressions for η(a, b) . . . . . . . . . . . . . . . . . . 193.2.3 Subproblems . . . . . . . . . . . . . . . . . . . . . . . 203.3 Exponential subproblems . . . . . . . . . . . . . . . . . . . . 22vTABLE OF CONTENTS3.4 Optimal Solutions [Loc16] . . . . . . . . . . . . . . . . . . . . 253.4.1 Solving the subproblems . . . . . . . . . . . . . . . . . 253.4.2 Functional forms . . . . . . . . . . . . . . . . . . . . . 263.4.3 Solutions . . . . . . . . . . . . . . . . . . . . . . . . . 263.5 Convex envelope: the maximum of all subproblems . . . . . . 293.6 Algorithmic Design . . . . . . . . . . . . . . . . . . . . . . . . 303.6.1 Input and Output Data Structures . . . . . . . . . . . 303.6.2 Main Algorithm . . . . . . . . . . . . . . . . . . . . . 333.6.3 Solving the subproblems . . . . . . . . . . . . . . . . . 373.7 Map Overlay problem . . . . . . . . . . . . . . . . . . . . . . 403.7.1 Data structures for polyhedral region . . . . . . . . . . 413.7.2 Sorting vertices in clockwise direction . . . . . . . . . 423.7.3 Map Overlay Algorithm . . . . . . . . . . . . . . . . . 45Chapter 4: Conjugate of a class of Convex Bivariate RationalFunctions over a polytope . . . . . . . . . . . . . . 494.1 Subdifferentials in the interior of the polytope . . . . . . . . . 504.2 Subdifferentials at the Vertices . . . . . . . . . . . . . . . . . 554.3 Subdifferentials on the edges . . . . . . . . . . . . . . . . . . 554.4 Structure of the conjugate domain . . . . . . . . . . . . . . . 594.5 Conjugate Expressions . . . . . . . . . . . . . . . . . . . . . . 654.5.1 Fractional forms . . . . . . . . . . . . . . . . . . . . . 65Chapter 5: Algorithmic computation of the Conjugate for aclass of Bivariate Rational functions over a poly-tope . . . . . . . . . . . . . . . . . . . . . . . . . . . . 735.1 Data Structures . . . . . . . . . . . . . . . . . . . . . . . . . . 735.1.1 Parabolic region . . . . . . . . . . . . . . . . . . . . . 745.1.2 Output Data Structure . . . . . . . . . . . . . . . . . 765.2 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 795.2.1 Main Algorithm . . . . . . . . . . . . . . . . . . . . . 805.2.2 Algorithm for edges . . . . . . . . . . . . . . . . . . . 815.3 Example 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 835.4 Example 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89Chapter 6: Conclusions and Future work . . . . . . . . . . . . 966.1 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 976.1.1 Maximum of all conjugates . . . . . . . . . . . . . . . 976.1.2 Conjugate of bivariate nonconvex PLQ Functions . . . 986.1.3 Convex envelope of bivariate nonconvex PLQ functions 98viTABLE OF CONTENTSBibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105Appendix A: The η(a, b) expressions for the bilinear functions . . . 106Appendix B: Solving the subproblems . . . . . . . . . . . . . . . . 108B.1 Case 1: Quadratic-Quadratic . . . . . . . . . . . . . . . . . . 108B.1.1 Solving the equality constraint . . . . . . . . . . . . . 109B.1.2 Solving the inequalities . . . . . . . . . . . . . . . . . 109B.1.3 Solutions . . . . . . . . . . . . . . . . . . . . . . . . . 111B.2 Case 2: Quadratic-Linear . . . . . . . . . . . . . . . . . . . . 112B.2.1 Solving the equality constraint . . . . . . . . . . . . . 113B.2.2 Solving the inequalities . . . . . . . . . . . . . . . . . 113B.2.3 Solutions . . . . . . . . . . . . . . . . . . . . . . . . . 115B.3 Case 3: Linear-Linear . . . . . . . . . . . . . . . . . . . . . . 116B.3.1 Solving the equality constraint . . . . . . . . . . . . . 116B.3.2 Solving the inequalities . . . . . . . . . . . . . . . . . 117B.3.3 Solutions . . . . . . . . . . . . . . . . . . . . . . . . . 118Appendix C: Conjugate expressions for a rational function over anedge . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119viiList of TablesTable 3.1 Different cases depending upon the eigenvalues . . . . 14Table 5.1 Conjugate domain for a rational function . . . . . . . 76Table 5.2 Conjugate domain for quadratic functions . . . . . . . 80Table 6.1 Observed intersections for computing the maximum ofall conjugates . . . . . . . . . . . . . . . . . . . . . . . 97viiiList of FiguresFigure 1.1 A univariate convex PLQ Function . . . . . . . . . . 2Figure 1.2 A non convex univariate PLQ Function . . . . . . . . 3Figure 1.3 f(x) =x22, its conjugate fâ(s) =s22, and the map-ping between the two . . . . . . . . . . . . . . . . . . 5Figure 1.4 Closed convex envelope (shown in Blue) of each pieceof a univariate nonconvex PLQ function . . . . . . . . 6Figure 2.1 A Polyhedral subdivision of R2 . . . . . . . . . . . . . 10Figure 3.1 Î·1âs polyhedral domain . . . . . . . . . . . . . . . . . 22Figure 3.2 η2(a, b)âs polyhedral domain . . . . . . . . . . . . . . 23Figure 3.3 Nine Subregions formed by overlaying domains of η1and η2 . . . . . . . . . . . . . . . . . . . . . . . . . . 23Figure 3.4 19 subregions, shown by different colors, formed byintersection of domains of 3 η(a, b)s belonging to aconvex edge. . . . . . . . . . . . . . . . . . . . . . . . 24Figure 3.5 Input Data Structure . . . . . . . . . . . . . . . . . . 31Figure 3.6 Output Data Structure . . . . . . . . . . . . . . . . . 32Figure 3.7 Space division by line y = mx+ c . . . . . . . . . . . 41Figure 3.8 Vertex Enumeration Problem . . . . . . . . . . . . . . 48Figure 4.1 Primal and dual mapping for quadratic function inExample 4.4 over lines parallel to x1 + x2 = 0. . . . 52Figure 4.2 Normal cone at a vertex v â P shown as the arrowsin the intersection of Red and Green region . . . . . . 56Figure 4.3 Polyhedral subdivision for Example 4.14 . . . . . . . 62Figure 4.4 Parabolic subdivision for r and P in Example 4.22 . . 64Figure 4.5 Parabolic subdivision for r and P from Example 4.5 . 66Figure 5.1 A subdifferential region . . . . . . . . . . . . . . . . . 74Figure 5.2 Region with only active constraints . . . . . . . . . . 75ixLIST OF FIGURESFigure 5.3 Proposed subdivision for storing a subdifferential region 75Figure 5.4 Output Data structure . . . . . . . . . . . . . . . . . 77Figure 5.5 Normal cone division lines for Example 5.3 . . . . . . 83Figure 5.6 Conjugate for Example 5.3 . . . . . . . . . . . . . . . 88Figure 5.7 Normal cone division lines for Example 5.4 . . . . . . 90Figure 5.8 Conjugate for Example 5.4 . . . . . . . . . . . . . . . 95xList of Notations andAbbreviationsR Real numbersÏ Empty set|A| Cardinality of set Adom(f) Effective domain of function fcl(g) Closed envelope of function gconvfP Convex envelope of f over a region Pconv Closed convex EnvelopeSnÃn Symmetric Matrix of size nà nint(P ) Interior of the set Pri(P ) Relative interior of set PPLQ Piecewise Linear QuadraticLFT Legendre-Fenchel transformfâ Conjugate of function fsupS The supremum of set Sinf S The infimum of set SNP (x) Normal cone of set P at xâf(x) Subdifferential of function f at xâf(x) Gradient of function f at xIP Indicator function of set PxiAcknowledgementsI would take this opportunity to acknowledge and express my sinceregratitude towards my supervisor Dr. Yves Lucet, for his constant guidance,encouragement, and precious advice throughout my graduate studies. I amvery grateful for his endless support, especially for his expertise, patienceand insightful suggestions during the course of my research. His attentionto details motivated me to pay closer attention to everything I work on, andhis positive outlook and confidence in my research inspired me and gave methe right confidence to dig deeper into my work.I would also like to extend my gratitude towards the members of my su-pervisory committee: Dr. Heinz Bauschke and Dr. Warren Hare for spendingtheir valuable time reading my thesis and providing constructive and insight-ful feedback. I am also thankful to all of my committee members for makingmy defense a wonderful moment of my life.I truly appreciate the support of my family and all of my friends andcolleagues of the CCA research group for all they meant to me during thecompletion of my research. Without them, this journey and the lab wouldnot have been such an enjoyable experience.I am also thankful to the University of British Columbia, Okanagan andthe Natural Sciences and Engineering Research Council of Canada (NSERC)for providing the financial assistance, without which this research would haveseemed like a distant dream.xiiDedicationTo my family,for their constant love and support.xiiiChapter 1IntroductionComputational convex analysis (CCA) focuses on creating efficient al-gorithms to compute fundamental transforms arising in the field of convexanalysis. Computing the convex envelope or biconjugate is the core oper-ation that bridges the domain of nonconvex analysis with convex analysis.The early idea of the computation of convex transforms can be traced backto [Mor65]. However, development of most of the algorithms in CCA beganwith the Fast Legendre Transform (FLT) in [Bre89] and was further devel-oped in [Cor96, Luc96] (and independently in [NV94, SAF92]). Later, thealgorithm was improved to the optimal linear worst-case time complexityin [Luc97]. Recent algorithms to compute the conjugate numerically havebeen based on either parameterization [HUL07], manipulation of graphs(GPH model) [GL11], or the computation of the Moreau envelope [Luc06].More complex operators such as the proximal average operator [BGLW08,BLT08] can be built by using a combination of addition, scalar multiplica-tion, and conjugacy operations. This has further motivated more researchin the field of computational convex analysis [LLT17, JKL11, LBT09].Computational convex analysis has numerous applications in a wide va-riety of fields, ranging from image processing and partial differential equa-tions to thermodynamics. An extensive study has been carried out in thesurvey [Luc10].1.1 PLQ FunctionsPiecewise Linear Quadratic (PLQ) functions play a significant role inthe approximation of functions. Compared to piecewise linear functions,PLQ functions capture the curve of the original function and are continu-ous but may not be smooth, however, they are second-order differentiablein the interior of each piece. PLQ functions have been well-known in thefield of convex analysis (the set of PLQ functions is closed under commonconvex operators [RW98]). Their applications in nontraditional statisticalestimations problems have been well studied in [CCHP17].11.2. Convex EnvelopeFigure 1.1: A univariate convex PLQ FunctionAn example of a convex univariate PLQ function is shown in Figure 1.1.A PLQ framework along with the existence of linear time algorithms forvarious convex transforms has been proposed in [LBT09, BLT08]. A detaileddescription of all the important transforms is presented in the study [Luc10].Computing the full graph of the convex hull of univariate nonconvex PLQfunctions using the biconjugate operation in optimal linear worst-case timecomplexity has been proposed in [GL10]. However, the area of bivariatePLQ functions, especially nonconvex, has not been explored much.1.2 Convex EnvelopeFor a function f defined over a region X, the process of computingthe pointwise supremum of all convex underestimators of f for region X iscalled the convex envelope and is denoted by convfX . For example, for theunivariate nonconvex PLQ functionf(x) =ï£±ï£´ï£²ï£´ï£³x2 x â¤ 01â (xâ 1)2 0 < x < 1x2 x â¥ 1shown in Figure 1.2, its convex envelope is given by21.2. Convex EnvelopeFigure 1.2: A non convex univariate PLQ Functionconvf =ï£±ï£´ï£²ï£´ï£³x2 x â¤ 0x 0 < x < 1x2 x â¥ 1and is shown earlier in Figure 1.1.Computing the convex envelope of a function is one of the fundamentaltask in mathematical optimization. Every nonconvex function shares thesame minima as its convex envelope and computing the minima of a convexfunction is much easier than computing the same for a nonconvex function.The task of computing the convex envelope of a function is very usefulthough not simple (computation of the convex envelope of a multilinearfunction over a unit hypercube has been proved to be NP-Hard in [Cra89]).However, the convex envelope of functions defined over a polytope P andrestricted by the vertices of P can be computed in finite time using a linearprogram [Tar04, Tar08]. A method to reduce the computation of convexenvelope of functions that are one lower dimension(Rnâ1) convex and haveindefinite Hessian to optimization problems in lower dimensions is discussedin [JMW08].Results for specific functions exist in the literature, a particular one isthe bilinear function. Any general bivariate nonconvex quadratic functioncan be linearly transformed to the sum of bilinear and a linear function.Some studies involving bilinear functions are31.3. Fenchel conjugateâ convex envelopes for bilinear functions over rectangles have been dis-cussed in [McC76] and validated in [AKF83];â the convex envelope over special polytopes (not containing edges withfinite positive slope) has been derived in [SA90];â the convex envelope of a bilinear function over a triangle containingexactly one edge with finite positive slope is presented in [Lin05];â the convex envelope over general triangles and triangulation of thepolytopes through doubly nonnegative matrices (both semidefinite andnonnegative) is presented in [AB10, Ans12].In [Loc14], it is shown that the analytical form of the convex envelope ofsome bivariate functions defined over polytopes can be computed by solvinga continuously differentiable convex problem and then later is extended toany bivariate function satisfying some assumption (see Assumption [Loc16]).They go on to prove that the convex envelope of any bivariate functionover a polytope satisfying the assumption is characterized by a polyhedralsubdivision, however, the form of the envelope is determined implicitly.1.3 Fenchel conjugateThe Fenchel conjugate of a function f : Rn â R âª {+â} is given byfâ(s) = supxâRn[ãs, xã â f(x)]where ãs, xã = sTx is the standard dot product. It is also known as theLegendre-Fenchel Transform or convex conjugate or simply conjugate. Itplays a significant role in duality and computing it is a key step in solvingthe dual optimization problem [Her16]. For a convex function f(x) =x22,the conjugate is the function itself in the dual domain. The mapping betweenthe primal and dual domain for this function is shown in Figure 1.3. Fora smooth convex univariate function, the slope of the affine underestimatorat a given point in the primal domain maps to a point in the dual domain,while the intercept multiplied by â1 becomes the value of the conjugate atthat slope and vice-versa.A method to compute the conjugate known as the fast Legendre trans-form was first introduced in [Bre89] and its applications have been studiedin [Cor96, Luc96]. A linear time algorithm was later introduced by Lucet tocompute the discrete Legendre transform [Luc97].41.4. Motivationy = 1*x-0.5r = 1*s -0.5Primal domain (x,y) Dual domain (s,r)1.510.5 0 -0.5 -1.5 -1 -0.5 0 0.5 10.5 0 -0.5-1-1.5 -1 -0.5 0 0.5 1 1.51Figure 1.3: f(x) =x22, its conjugate fâ(s) =s22, and the mapping betweenthe twoComputation of the conjugate of convex univariate PLQ functions havebeen well studied in the literature and linear time algorithms have beendeveloped in [GL13, GJL14] both for the PLQ and the GPH models. Re-cently, a linear time algorithm to compute the conjugate of convex bivariatePLQ functions using an entity graph (storing the entities of the domain asa graph) has been proposed in [HL18].1.4 MotivationLet f : Rn â R âª {+â} be a piecewise function, that is,f(x) =ï£±ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£³f1(x) x â P1f2(x) x â P2· · · · · Â·fN (x) x â PN .From [HUL93, Theorem 2.4.1], we have(infifi)â = supifâi . (1.1)51.4. MotivationFigure 1.4: Closed convex envelope (shown in Blue) of each piece of a uni-variate nonconvex PLQ functionAlso, from [HUL93, Proposition 2.6.1],conv(infi(fi + IPi)) = conv( infi[conv(fi + IPi)] ) (1.2)where IPi is the indicator function for Pi. Soconv(infi(fi + IPi)) = conv(infi[conv(fi + IPi)]) using (1.2)= ((infi[conv(fi + IPi)])â)â= ( supi[conv(fi + IPi)]â )â using (1.1).Hence the closed convex envelope of a PLQ function can be computedin four steps:Step 1 Compute the convex envelope of each pieceStep 2 Compute the conjugate of the convex envelope of each piece61.4. MotivationStep 3 Compute the maximum of all the conjugatesStep 4 Compute the conjugate of the function obtained in Step 3 to get thebiconjugate (or the closed convex envelope) of the original function.Note that computing the closed convex envelope of each piece does not resultin the closed convex envelope of the whole function. This is illustrated inFigure 1.4.In this thesis, we begin with computing the convex envelope of eachpiece using the method proposed in [Loc16] in Chapter 3. After that, wecompute the conjugate of each piece in Chapter 4. Chapter 5 covers thealgorithmic design and examples for computing the conjugate of a convexrational function over a polytope. Step 3 and 4 are discussed in Chapter 6as future work.7Chapter 2Preliminaries and NotationsThis section covers the basic background to facilitate proper understand-ing of the rest of the thesis.Definition 2.1. Convex Set. A set S is convex if âx, y â S and λ â [0, 1],λx+ (1â λ)y â S.Definition 2.2. Cardinality of a finite set. Cardinality of a finite set is thenumber of elements in the set.For example, let A = {1, 5, 7}, then its cardinality |A|= 3.Definition 2.3. Convex Function. A function f : Rn â Râª{+â} is convexif and only if âx, y â dom(f) = {x : f(x) â R} and λ â (0, 1),f(λx+ (1â λ)y) â¤ λf(x) + (1â λ)f(y).It is strictly convex whenf(λx+ (1â λ)y) < λf(x) + (1â λ)f(y)for x, y â dom(f), x 6= y and λ â (0, 1).Definition 2.4. Lower Semi-continuous/Closed function [HUL13, Defini-tion 3.2.1]. A function f : Rn â R âª {+â} is said to be closed iflim infxâx0f(x) â¥ f(x0) for all x0 â Rn.For any non closed function g : Rn â Râª {ââ,+â}, we note cl(g) theclosed envelope of g, that is, the largest closed function that is lower thang. We denote the set of closed convex functions on Rn as convRn. All thefunctions involved here on are in convRn.Definition 2.5. Convex Envelope. The convex envelope of a function fover a region P , noted conv fP , is the pointwise supremum of all the convexunderestimators of f over P . It is defined asconvfP (x, y) = sup cf(xâ², yâ²)â [a(xâ² â x) + b(yâ² â y) + c] â¥ 0, â(xâ², yâ²) â P8Chapter 2. Preliminaries and NotationsThe optimal value (aâ, bâ, câ) defines the supporting hyperplane for thefunction f at a point (x, y) â P .Definition 2.6. Symmetric Matrix. A symmetric matrix, A, is a squarematrix which is equal to its transpose, that is, A = AT . The set of symmetricmatrices is denoted by S.Definition 2.7. Positive Semidefinite Matrix. An nà n symmetric matrixM is called positive semidefinite if for any x 6= 0 â Rn, xTMx â¥ 0. WhenxTMx > 0, then M is a positive definite matrix.Definition 2.8. Indefinite Matrix. A symmetric matrix, M , is called in-definite if it contains at least two nonzero eigenvalues of opposite signs.For a function f â C1, an indefinite Hessian signifies f is nonconvex.Definition 2.9. Halfspace. A halfspace is an unbounded region formed bythe set of points satisfying a linear inequality. Mathematically,Hs = {x â Rn : aTx â¤ b}where a â Rn, a 6= 0 and b â R.Definition 2.10. Hyperplane(Affine). A hyperplane is formed by the set ofpoints statisfying a linear equality. Mathematically,Hp = {x â Rn : aTx = b}where a â Rn, a 6= 0, b â R.Definition 2.11. Polytope. A polytope P is a bounded region formed bythe intersection of a finite number of halfspaces. It can be written asP = {x â Rn : Ax â¤ b},where A â RmÃn, b â Rm and m is the number of halfspaces. A polytope isalways convex.A polytope in R2 can be divided into three entities namely vertices, edges,and the interior. A vertex is formed by the intersection of edges. An edge isa line segment on the boundary of the polytope joining two vertices, whilethe interior of a polytope P can be written as int(P ) = {x â Rn : Ax < b}.9Chapter 2. Preliminaries and NotationsFigure 2.1: A Polyhedral subdivision of R2Definition 2.12. Polyhedral Subdivision. A convex set R =âiâ{1,...,m}Ri,R â R2, defined as the union of a finite number of polyhedral regions is saidto have a polyhedral subdivision if for any j, k â {1, · · · ,m}, j 6= k,RjâRkis either Ï or is contained in a line.Figure 2.1 shows an example of polyhedral subdivision of R2.Definition 2.13. Piecewise Linear-Quadratic function. A Piecewise Linear-Quadratic(PLQ) function is characterized by a polyhedral subdivision suchthat over each member of the subdivision, it has either linear or quadraticform. Mathematically, a PLQ functionf(x) =ï£±ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£³12xTA1x+ b1x+ δ1 x â P112xTA2x+ b2x+ δ2 x â P2· · · · · Â·12xTANx+ bNx+ δN x â PN .where x â Rn, Ai â SnÃn, bi â Rn, δi â R and Pi are polyhedral sets.Definition 2.14. Legendre-Fenchel transform. For a function f : Rn âR âª {+â}, its Legendre-Fenchel transform(LFT), fâ, is defined as10Chapter 2. Preliminaries and Notationsfâ(s) = supxâRn[ãs, xã â f(x)]where ãs, xã = sTx is the standard dot product.Definition 2.15. Biconjugate. For a function f : Rn â R âª {+â}, itsbiconjugate is given as (fâ)â.When f is convex (fâ)â = cl(f), otherwise (fâ)â is the closed convexenvelope of f .Definition 2.16. Cone. A set C is a cone if is contains 0 and for everyx â C, λx â C for all λ â¥ 0.Definition 2.17. Relative Interior [HUL13, Definition 2.1.1]. The relativeinterior of a convex set C â Rn, denoted by ri(C), is defined as the interiorwithin the affine hull of C. Mathematically, x â ri(C) if and only if x âaff(C) and âÎ´ > 0 such that aff(C) â©B(x, δ) â C where aff(C) is the affinehull of C and B(x, δ) is a ball of radius δ centered on x.Definition 2.18. Normal Cone. For a set P , the normal Cone at a pointx â P is given byNP (x) = {s : ãs, y â xã â¤ 0, for all y â P}.For all x â int(P ), NP (x) = {0}. When P is a 2D polytope, for any edgeE = {x : L(x) â¤ 0, l1 â¤ x1 â¤ u1} with L(x) = ax+ by+ c between vertices land u, the normal cone at any x â ri(E), NP (x) = {s : s = Î»âL(x), λ â¥ 0}.Definition 2.19. Subdifferential of a function. The subdifferential of afunction f at any x â dom(f), given by âf(x) isâf(x) = {s : f(y) â¥ f(x) + ãs, y â xã, for all y â dom(f)}.Every element of âf(x) is called a subgradient of f at x. When f is differ-entiable at x then âf(x) = {âf(x)}.Definition 2.20. Parabola (or Parabolic curve). A parabola is a two dimen-sional planar curve whose points (x, y) satisfy the equation ax2+bxy+cy2+dx+ ey+f = 0 where b2â4ac = 0. A line is a parabola with a = b = c = 0.11Chapter 2. Preliminaries and NotationsDefinition 2.21. Parabolic region. A parabolic region, Pr â R2, is formedby the intersection of a finite number of parabolic inequalities. It can bedefined asPr = {x â R2 : aix21 + bix1x2 + cix22 + dix1 + eix2 + fi â¤ 0, i â {1, · · · , k},b2i â 4aici = 0}.For Cp(x) = ax21 + bx1x2 + cx22 + dx1 + ex2 + f with b2 â 4ac = 0, the setPr = {x â R2 : Cp(x) â¤ 0} is convex, but Ps = {x â R2 : Cp(x) â¥ 0} is not.Definition 2.22. Parabolic subdivision. A convex set R =âiâ{1,...,m}Ri,R â R2, defined as the union of a finite number of parabolic regions is saidto have a parabolic subdivision if for any j, k â {1, · · · ,m}, j 6= k,RjâRkis either Ï or is contained in a parabola.12Chapter 3Convex Envelope ofBivariate QuadraticFunctions over a polytopeComputing the closed convex envelope of a 2D PLQ function requirescomputation of the closed convex envelope of each quadratic piece. For thispurpose, we utilize the method proposed in [Loc16]. We start with problemformulation, reduce it to subproblems, and then obtain the optimal solutionfor each. The closed convex envelope then comes out to be the maximum ofall such solutions over the given polytope.In this chapter, we formulate the mathematical form of the closed convexenvelope of each piece of a bivariate nonconvex PLQ function using [Loc16].After that, we introduce the data structures and algorithms involving various2D geometry manipulations to compute the closed convex envelope of abivariate nonconvex quadratic function.Before we present the method proposed in [Loc16] and our algorithmicapproach, let us discuss a technique to reduce the general expression of anynonconvex 2D quadratic function to a sum of bilinear and affine functions.3.1 Quadratic ReductionUsing linear algebra and eigenvalue decomposition, we can reduce anyquadratic function to a sum of a bilinear function and an affine function.This reduction is required to simplify our process of solving the subprob-lems, which, otherwise, would have led to solving quartic equations in mostof the cases. The other benefit of this reduction is that solving the subprob-lems only for the bilinear term, as compared to the general quadratic form,reduces the size, complexity, and number of special cases to be dealt withwhen creating an algorithmic solution for the problem.A general quadratic expression in R2 can be written asf(x) = xTAx+ bTx+ δ (3.1)133.1. Quadratic ReductionTable 3.1: Different cases depending upon the eigenvaluesEigenvalues Quadratic function type Convex envelope overa polytopeλ1 â¥ 0, λ2 â¥ 0 convex Function itselfλ1 > 0, λ2 > 0 strictly convex Function itselfλ1 â¤ 0, λ2 â¤ 0 concave Piecewise Linearλ1 < 0, λ2 < 0 strictly concave Piecewise Linearλ1 > 0, λ2 < 0 nonconvex and nonconcave Computed using [Loc16]where A â S2Ã2, b â R2, δ â R and S is the set of all symmetric matrices(see Definition 2.6). When f is neither convex or concave, A is an indefinitematrix.Let f1(x) = xTAx and f2(x) = bTx+ δ. From [HUL13, p. 71], we haveconvf1 + convf2 â¤ conv(f1 + f2),and because f2 is affine,convf1 + convf2 = conv(f1 + f2).An affine function is already convex, soconvf = convf1 + bTx+ δ. (3.2)For f1(x) = xTAx, by using eigenvalue decompositionA = UT[λ1 00 λ2]U, (3.3)where λ1 and λ2 are the eigenvalues, and U is the matrix of eigenvectors ofA.Let xË = Ux, thenf1(Uâ1xË) = g1(xË) = xËT[Î»1 00 Î»2]xËorg1(xË) = Î»1xËT[1 00 âÎ±]xË,(3.4)with α = âÎ»2/Î»1. When f(x) is convex, α â¤ 0. Since we are interestedin the case when f(x) is neither convex nor nonconcave, λ1 and λ2 would143.1. Quadratic Reductioncarry different signs, hence α > 0. Table 3.1 enumerates different possiblecases based on different eigenvalues.Now [1 00 âÎ±]=[1 1âÎ± ââÎ±] [0 11 0] [1 1âÎ± ââÎ±]T. (3.5)Let x¯ =[1 1âÎ± ââÎ±]TxË, so by using (3.5) in (3.4), we geth1(x¯) = λ1x¯T[0 11 0]x¯.Hence for x¯ = [x¯1 x¯2]T ,h1(x¯) = λ1x¯1x¯2 (3.6)which is our bilinear term.All the linear transformations performed in (3.3)-(3.6) are invertible asU is invertible and matrix[1 1âÎ± ââÎ±]is clearly invertible with inverse12âÎ±[âÎ± 1âÎ± â1].Proposition 3.1. For a bivariate quadratic nonlinear function f ofform (3.1), a polytope P , let f(x) = f1(x) + f2(x) with f1(x) = xTAx andf2(x) = bTx+ δ, then for λ1 > 0 and h1(x) = x1x2, conv fP = λ1conv h1P â¦Mâ1 + conv f2P where M is a 2Ã 2 invertible matrix.Proof. From [HUL13, p. 71], for any function f and λ > 0convλf = λ convf.So for (3.6), convf1 = λ1 convh1. The value λ1 > 0 can be made sure byalways assigning the positive eigenvalue to λ1.Now, for M =[1 1âÎ± ââÎ±]TU , h1(x¯) = f1(M â¦ x). Soconvf1 = convh1 â¦Mâ1. (3.7)Finally, (3.2) becomesconvf = λ1convh1 â¦Mâ1 + bTx+ δ.153.2. Problem Formulation[Loc16]Hence the closed convex envelope of any nonconvex quadratic functionf defined over a polytope P can be obtained byconvfP = λ1convh1P â¦Mâ1 + convf2P .Example 3.2. For the quadratic function f(x) = x[3 â1â1 â5]xT+[2,â2]Tx+4, its quadratic term gets reduced to h1(x¯) â 3.123 · x¯1x¯2 with transfor-mation matrix M â[â0.8360 â1.14901.3934 â1.1490].Since the problem of computing the closed convex envelope of a quadraticfunction reduces to computing the closed convex envelope of a bilinear func-tion, from now on we focus on formulating the problems, its solutions, andthe algorithmic design for bilinear function only. Also, for consistency of thischapter with [Loc16], we use variables x, y â R instead of x = [x1, x2] â R2in the algebraic expressions.3.2 Problem Formulation[Loc16]In this section, we will reformulate the problem of computing the closedconvex envelope into a group of subproblems as presented in [Loc16]. Themaximum of the solutions of such subproblems, when computed over thegiven polytope, leads to a polyhedral subdivision of the polytope such thatover each such division, the closed convex envelope belongs to a particularclass of bivariate rational functions.3.2.1 Problem FormulationThis section covers a new formulation of the closed convex envelopeproblem presented before in Definition 2.5. The closed convex envelope ofany function f at a point (x0, y0) â P , where P â R2 is a convex set, isgiven by an optimization problem in a, b, c with infinite number of linearconstraints asconvfP (x0, y0) = max cf(x, y)â [a(xâ x0) + b(y â y0) + c] â¥ 0 â(x, y) â P.The optimal solution (aâ, bâ, câ) of this problem defines the supporting hy-perplane to the closed convex envelope of f at point (x0, y0) over P .163.2. Problem Formulation[Loc16]The infinite number of constraints can be replaced with a constraint asconvfP (x0, y0) = max cmin(x,y)âPf(x, y)â [a(xâ x0) + b(y â y0)] â¥ c.Let us introduce the following assumption,Assumption 3.2.1. For a polytope P â R2 and function f â C2,â the Hessian of f is indefinite in the interior of P , andâ along each edge of P , f is either strictly convex or concave.This first part of the assumption is sufficiently reasonable as we want tocompute the closed convex envelope of a function which is neither convexnor concave. However the second condition is restrictive and requires carefulselection of the polytope P , but it plays a crucial role in the new formula-tion of the convex envelope problem. Most of the 2D nonconvex nonlinearquadratic functions over any polytope P satisfy these assumptions. For thequadratic functions that do not statisfy these assumptions, a change in thepolytope can be used to make them valid under our assumptions.Now under the above assumptions, the minimum of f(x, y)â [a(xâx0)+b(yâ y0)] â¥ c will always be attained either at a vertex of P or over an edgesuch that the constriction of f along the edge is a strictly convex function.So the problem of computing the closed convex envelope of f over a polytopeP can then be written asconvfP (x, y) = max cf(xvi , yvi)â [a(xvi â x) + b(yvi â y)] â¥ c, â i â V â²(P )min(xâ²,yâ²)â[vj1 ,vj2 ]f(xâ², yâ²)â [a(xâ² â x) + b(yâ² â y))] â¥ c,âj âEâ²(P )(3.8)whereâ Eâ²(P ) is the index set of edges where f is strictly convex,â âj â Eâ²(P ), vertices vj1 and vj2 are the end points of the edge,â V â²(P ) is the set of all the remaining vertex indices not belonging inEâ²(P ),173.2. Problem Formulation[Loc16]â âi â V â²(P ), (xvi , yvi) are the coordinates of the corresponding vertex.Example 3.3. Let f(x, y) = xy and P be a polytope with verticesv1 = (0, 0), v2 = (1, 2), v3 = (2, 2), v4 = (2, 1).The equation of lines constituting the boundaries of polytope are e1 := y =2x from v1 to v2, e2 := y = 2 from v2 to v3, e3 := x = 2 from v3 to v4, ande4 := y = x/2 from v4 to v1 respectively. Since f is strictly convex alongthe lines y = 2x and y = x/2 only, so Eâ²(P ) = {1, 4} and V â²(P ) = {3}.A detailed formulation of the above steps has been presented in [LS14,Sec. 2]. In the above formulation, each constraint either belongs to an edgeor a vertex of the polytope P . However, the constraint belonging to an edgestill remains a minimization problem in itself.For j â Eâ²(P ), a convex edge with equation of line as yâ² = mjx + qj ,let fj(xâ²) = f(xâ²,mjx+ qj) be the restriction of f along the edge. Now theminimization problem related to this edge can be written asminxâ²â[xvj1 ,xvj2 ]fj(xâ²)â [a(xâ² â x) + b(mjxâ² + qj â y)] â¥ c. (3.9)Let sj(a, b) = arg min{(3.9)}. If the function is decreasing (resp. in-creasing), sj(a, b) = +â(resp. ââ), otherwisef â²j(sj(a, b)) = a+ bmj . (3.10)Since fj is strictly convex along [xvj1 , xvj2 ], we have fâ²j(xvj2 ) â¥ f â²j(xvj1 ).Now the minima will either lie in the interval or on the boundaries. Let usdefine the three regionsDâj = {(a, b) : a+ bmj â¤ f â²j(xvj1 )}Dj = {(a, b) : f â²j(xvj1 ) â¤ a+ bmj â¤ f â²j(xvj2 )}D+j = {(a, b) : a+ bmj â¥ f â²j(xvj2 )}.(3.11)The arg min{(3.9)} isxj(a, b) =ï£±ï£´ï£²ï£´ï£³xvj1 if (a, b) â Dâjsj(a, b) if (a, b) â Djxvj2 if (a, b) â D+j .(3.12)Correspondingly,yj(a, b) = mjxj(a, b) + qj .183.2. Problem Formulation[Loc16]Also for the constraints belonging to the remaining vertex indices i âV â²(P ),xi(a, b) = xvi and yi(a, b) = yvi .So âk â V â²(P ) âª Eâ²(P ), we defineηk(a, b) = f(xk(a, b), yk(a, b))â axk(a, b)â byk(a, b). (3.13)By using (3.13), the closed convex envelope problem (3.8) becomesconvfP (x, y) = max cηk(a, b) + ax+ by â¥ c,â k â V â²(P ) âª Eâ²(P ).(3.14)The problem (3.14) is a simplified version of (3.8) where ηk can taketwo different forms depending upon if it is a convex edge or a vertex of thepolytope.3.2.2 Expressions for η(a, b)The bilinear function, f(x, y) = xy, is strictly convex on an edge y =mjx + qj if mj > 0 (as fâ²â²j (x) = 2mj and to be strictly convex fâ²â²j (x) â¥ 0).From (3.12), over such edgesj(a, b) =a+ bmj â qj2mj.Soâj â Eâ²(P ),ηj(a, b) =ï£±ï£²ï£³â(a+mjbâ qj)24mjâ bqj (a, b) â Djâaxj(a, b)â byj(a, b) + xj(a, b)yj(a, b) (a, b) â {Dâj , D+j },and, âi â V â²(P ) and vertex (xvi , yvi),ηi(a, b) = âaxvi â byvi + xviyvi .(3.15)The derivation of the above expressions is presented in Appendix A.Example 3.4. From Example 3.3, we have two convex edges, e1 and e4with lines y = 2x and y = x/2 between vertices v1 = (0, 0) to v2 = (1, 2)193.2. Problem Formulation[Loc16]and v1 = (0, 0) to v4 = (2, 1) respectively. Now for the edge e1,η1(a, b) =ï£±ï£´ï£´ï£²ï£´ï£´ï£³0 a+ 2b â¤ 0â(a+ 2b)280 â¤ a+ 2b â¤ 42â aâ 2b 4 â¤ a+ 2b.Similarly, for edge e4η2(a, b) =ï£±ï£´ï£´ï£²ï£´ï£´ï£³0 a+ b2 â¤ 0â(a+ b/2)220 â¤ a+ b2 â¤ 22â 2aâ b 2 â¤ a+ b2 .The indexing of η functions does not matter as our main objective is toobtain the constraint set only. Finally, for vertex v3η3(a, b) = 4â 2aâ 2b.3.2.3 SubproblemsFrom [Loc16, lemma 3.1], at any point in the polytope not lying on thestrictly convex edge or vertices, the solution of problem (3.14) has at leasttwo active constraints. So the optimization problem (3.14) can be reducedto subproblems for each pair of η functions. Hence, â(h,w) â Eâ²(P )âV â²(P )and h 6= w,convfP (x, y) = maxh,wmaxa,bηh(a, b) + ax+ byηh(a, b) = ηw(a, b)ηh(a, b) â¤ ηk(a, b), k â Eâ²(P ) âª V â²(P ) {h,w}(3.16)Since η is a piecewise defined function with a polyhedral subdivision,solving the subproblems for each pair is the same as solving them over eachmember of the resulting polyhedral subdivision obtained by overlaying allthe polyhedral subdivisions for all η functions.When arranged according to the subdivided region, the form of the sub-203.2. Problem Formulation[Loc16]problems becomemax ηh(a, b) + ax+ byηh(a, b) = ηw(a, b)ηh(a, b) â¤ ηk(a, b), k â Eâ²(P ) âª V â²(P ) {h,w}(a, b) â Sr(3.17)whereâ index r of subproblem ranges in interval 1, ..., 3|Eâ²(P )| where |Eâ²(P )| isthe cardinality of set Eâ²(P )(see Definition 2.2),â Each polyhedral set Sr is defined asSr = â©jâEâ²(P )Îjwhere Îj â {Dâj , Dj , D+j } i.e all possible combinations of Dâj , Dj , D+jand âj â Eâ²(P ).A set S = Sr is the collection of all the polyhedral set such that overeach Sr, there is an unique combination of pieces of η functions. Also,âiâ{1,...,|S|} Sr = R2.Example 3.5. From Example 3.4, the set S = {{a + b/2 â¤ 0, a + 2b â¤0}, {a+ b/2 â¤ 0, 0 â¤ a+ 2b â¤ 4}, {a+ b/2 â¤ 0, a+ 2b â¥ 4}, · · · , {a+ b/2 â¥2, a+ 2b â¥ 4}}.Now the subproblem defined over subregion S2 ismax 0 + ax+ by0 = â(a+ 2b)280 â¤ 4â 2aâ 2b(a, b) â S2 = {a+ b/2 â¤ 0, 0 â¤ a+ 2b â¤ 4}.Similarly, the subproblem defined over subregion S5 ismax â (a+ 2b)28+ ax+ byâ (a+ 2b)28= â(a+ b/2)22â (a+ 2b)28â¤ 4â 2aâ 2b(a, b) â S5 = {0 â¤ a+ b/2 â¤ 2, 0 â¤ a+ 2b â¤ 4}.213.3. Exponential subproblemsFigure 3.1: Î·1âs polyhedral domain3.3 Exponential subproblemsWe rearranged the subproblems from (3.16) to (3.17) by using the infor-mation that each η that belongs to a convex edge, is a piecewise functionwith exactly 3 pieces defined over a polyhedral set. Solving the subprob-lems for each pair (h,w) would eventually lead to solving them for eachintersecting polyhedral set Sr. Now, the number of subproblem ranges in{1, ..., 3|Eâ²(P )|} because with each η the aâ b space could potentially be di-vided into 3 times the number of current divisions. This is best illustratedwith the help of Figures 3.1- 3.3.Figures 3.1 and 3.2 show the polyhedral domains for some η1(a, b) andη2(a, b). For the bilinear function, the domain of ηj(a, b) âj â Eâ²(P ) willalways be divided by two parallel hyperplanes of form a+mjb = fâ²(xvl) forl â {1, 2}, where v1 and v2 are the end points of the respective convex edge,and mj > 0.Now solving the subproblems in (3.16) for η1(a, b) and η2(a, b) is the sameas solving all the subproblems defined over each polyhedral set in Figure 3.3.Since every polyhedral set contains exactly one piece from each η(a, b), theformulation (3.17) is the rearrangement of subproblems according to theirsubregions.These Subregions are formed as a result of the intersection between re-gions Dâj , D+j and Dj for all ηj . The number of such possible regions are3|Eâ²(P )|, however, out of all the 3|Eâ²(P )| combinations, few will always beinvalid due to the restriction mj > 0. For the new subdivision to have 3times the number of current members, each polyhedral set of the addingsubdivision needs to have a valid intersection with all the members of the223.3. Exponential subproblemsa + m2b â¤ fâ²2(xv1) a + m2b â¥ fâ²2(xv2)â¤ fâ²2(xv2)fâ²2(xv1) â¤ a + m2bFigure 3.2: η2(a, b)âs polyhedral domainFigure 3.3: Nine Subregions formed by overlaying domains of η1 and η2233.3. Exponential subproblemsFigure 3.4: 19 subregions, shown by different colors, formed by intersectionof domains of 3 η(a, b)s belonging to a convex edge.current subdivision. This only holds upto 2 layers as each new layer of sub-division is of polyhedral form, however, 3|Eâ²(P )| still serves as a good upperbound on the number of subregions. This exponential bound also leads toa O (3N) time complexity algorithm, as any algorithmic approach to thisproblem needs to generate all the subregions, which is growing exponentiallywith the number of convex edges. A result involving 19 subregions formedby intersecting 3 ηj(a, b) for j â Eâ²(P ) is shown in Figure 3.4.The problem of obtaining all the valid subregions belong to a special classof map overlay problems and an algorithm to obtain all such subregions isdiscussed in section 3.7. Since there is an exponential number of subregions,the time complexity of the algorithm is exponential as well.Although the number of subregions is exponential, the number of sub-problems are much more than that. The ηi(a, b) âi â V â²(P ) do not con-tribute to any of the division of subregions as they are defined over the fulldomain. However, they add up to the number of subproblems to be solved243.4. Optimal Solutions [Loc16]over a subregion. For N = |Eâ²(P )| and M = |V â²(P )|, we have to solvethe subproblems for each η pair which leads to(N+M2)problems over eachsubregion. Finally, the total subproblems are bounded above by(N+M2)3Nwhich increases the time complexity of the solving algorithm further thanexponential. The best strategy, in this regard, is to choose a polytope thatcontains a low number of convex edges.3.4 Optimal Solutions [Loc16]In this section, we introduce the method proposed in [Loc16, Section4], and then formulate the possible forms of the optimal solutions for thesubproblems obtained in Section 3.2.3.3.4.1 Solving the subproblemsThe method proposed can be summarized as followsStep 1 Use the equality constraint to find the relation between the variableswhich leads to eliminating one variable from the subproblem. Thenuse this relation to reduce the subregion to one dimension.Step 2 Using the relation in step 1, obtain the feasible region by solving allthe inequalities. The feasible region will either get reduced to a unionof finite number of one-dimensional intervals or become infeasible. Inthe former case, the subproblem gets further divided into subproblemsover each interval.Step 3 Each one-dimensional problem is then solved to obtain the closed formsolutions with a polyhedral subdivision of R2.Step 4 After obtaining the solutions for each subproblem, the closed convexenvelope is the maximum of all such optimal solutions over the poly-tope. This again, leads to a further subdivision of the polytope suchthat over each member of the subdivision, the closed convex envelopehas one of the explicit forms discussed in Section 3.4.2.For the bilinear function, η has either a quadratic or a linear form as de-fined in (3.15). Solving the equality constraint in Step 1 then leads to threecases: quadratic-quadratic(qq), quadratic-linear(ql), and linear-linear(ll).For qq and ll, the relation between a and b comes out to be linear of forma = α1b+ α0, for α1, α0 â R. The case qq leads to two such linear relations253.4. Optimal Solutions [Loc16]resulting from the roots of a quadratic function. The case ql is special andrequires a new variable Z such that b = β2Z2 + β1Z + β0 for βi â R, anda = Z âmjb. These relations are then used to reduce the subregions, Sr, toa union of one dimension intervals.Step 2 then leads to solving inequalities for all the remaining η func-tions over each one-dimensional region. Step 3 leads to the formulation andcomputing the solution of the one dimension problem, while Step 4 involvescomputing the maximum of all the solutions for one dimension optimizationproblems. A detailed formulation of steps 1-3 is presented in Appendix B.3.4.2 Functional formsLet us introduce the following functional formsfr(x, y) =ξ21(x, y)ξ2(x, y)+ ξ0(x, y), (3.18)fq(x, y) =ξ21(x, y)ζ+ ξ0(x, y), ζ > 0 (3.19)fl(x, y) = ζξ1(x, y) + ξ0(x, y), (3.20)where ξi(x, y) are linear functions in x and y, and ζ â R.The functional forms fq and fl are just special cases of fr and can be ob-tained by setting coefficients of x, y in ξ2(x, y) to zero, and the constant termto ζ and 1/ζ respectively. Let L = {(x, y) : α10x + α01y + α00 = 0} whereα01, α10, α00 â R. In the next section, we will observe that the solutions toall the subproblems have one of the three functional forms mentioned above.3.4.3 SolutionsAfter solving the equality and all the inequalities in Step 1 and 2, eachsubproblem gets reduced to the formmax ηh(a, b) + ax+ bydl â¤ (b or Z) â¤ du(3.21)By substituting the relation between a and b in (3.21), we getmax â ξ2(x, y)u2 + 2uξ1(x, y) + ξ0(x, y)dl â¤ u â¤ du(3.22)where u = b or Z (depending upon the case) and ξ2, ξ1 and ξ0 are linearfunctions in x and y. Since (3.22) is differentiable and concave in u, tocompute the maximum we compute the critical points.263.4. Optimal Solutions [Loc16]Case 1: ηh(a, b) quadratic, ηw(a, b) quadraticWhen both η(a, b) functions are quadratic, the equality constraint re-duces to the linear relationsa = α1b± α0.After substituting one of the above relation in (3.21), we obtainmax âÎ¾2b2 + 2bξ1(x, y) + ξ0(x, y)dl â¤ b â¤ du(3.23)where ξ1, ξ0 â L, ξ2 > 0, andξ2 =(α1 +mh)24mhâ¥ 0,ξ1(x, y) = â(Î±1 +mh)(α0 â qh)2mh+ y â qh + α1x,ξ0(x, y) = â(Î±0 â qh)24mh+ α0x.Since (3.23) is differentiable and strictly concave in b, its optimal valueis ï£±ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£³Î¾21(x, y)ξ2+ ξ0(x, y) if dlξ2 < ξ1(x, y) < duÎ¾2,âÎ¾2d2l + 2dlξ1(x, y) + ξ0(x, y) if dlξ2 â¥ ξ1(x, y),âÎ¾2d2u + 2duξ1(x, y) + ξ0(x, y) otherwise.(3.24)For quadratic-quadratic case, the potential convex envelopeâs closed formexpression is defined by a polyhedral subdivision with three pieces (3.24))such that over each piece, it has form (3.19) or (3.20).Case 2: ηh(a, b) quadratic, ηw(a, b) linearIn this case, we define a variable Z = a + mhb such that the equalityconstraint leads tob = β2Z2 + β1Z + β0a = Z âmhb.Now the optimization problem becomesmax â ξ2(x, y)Z2 + 2Zξ1(x, y) + ξ0(x, y)dl â¤ Z â¤ du(3.25)273.4. Optimal Solutions [Loc16]where ξ2, ξ1, ξ0 â L andξ2(x, y) = Î²2xâ β1y + qhβ2 + 14mh,ξ1(x, y) =12((1â β1)x+ β1y + qh2mhâ qhβ1),ξ0(x, y) = âÎ²0x+ β0y â q2h4mhâ qhβ0.Since, in this case, the concavity of the objective function depends uponξ2, when ξ2(x, y) > 0, the optimal value of the subproblem isï£±ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£³Î¾21(x, y)ξ2(x, y)+ ξ0(x, y) if dlξ2(x, y) < ξ1(x, y) < duξ2(x, y),âÎ¾2(x, y)d2l + 2dlξ1(x, y) + ξ0(x, y) if dlξ2(x, y) â¥ ξ1(x, y),âÎ¾2(x, y)d2u + 2duξ1(x, y) + ξ0(x, y) otherwise.(3.26)and when ξ2(x, y) â¤ 0,max{âÎ¾2(x, y)d2l +2dlξ1(x, y)+ξ0(x, y),âÎ¾2(x, y)d2u+2duξ1(x, y)+ξ0(x, y)}.(3.27)The expression (3.27) may further divide the region ξ2(x, y) â¤ 0 into twosubregions if the intersecting line of both functions passes through ξ2(x, y) â¤0.Finally, in this case, the potential closed convex envelope function ischaracterized by a polyhedral subdivision with 4(or 5) pieces such that overeach piece it has form (3.18) or (3.20).Case 3: ηh(a, b) linear, ηw(a, b) linearIn this case, the equality constraint leads to a linear relationa = α1b+ α0.Since both η(a, b) functions are linear, the subproblem reduces tomax 2bξ1(x, y) + ξ0(x, y)dl â¤ b â¤ du(3.28)where ξ1, ξ0 â L withξ1(x, y) = α1x+ y â α1xh â yh,ξ0(x, y) = α0x+ xhyh â α0xh.283.5. Convex envelope: the maximum of all subproblemsSince the objective function is linear, the optimal value lies at the bound-aries, that is,max {dlξ1(x, y) + ξ0(x, y), duξ1(x, y) + ξ0(x, y)}.The solution, in this case, leads to a polyhedral subdivision with onlytwo pieces such that over each piece, the potential convex envelope hasform (3.20).3.5 Convex envelope: the maximum of allsubproblemsThe solutions of all the subproblems can be generalized as a piecewisefunction of form (3.18)-(3.20) with a polyhedrally divided domain in R2.Now the closed convex envelope is the pointwise maximum of all suchsolutions, that is,convfP = maxSâSuS(x, y)where S is the set of all subproblems and uS is the solution of each.Every solution is defined over the full domain R2 and each solution pro-vides at least two candidate expressions for the closed convex envelope.Computing the closed convex envelope leads to solving a map overlay prob-lem with each solution as a different layer. Since we are only interested inthe functions defined over the polytope P , the contribution from each sub-problem reduces to at least one. Hence for N number of valid subproblems,we need to find the maximum of N expressions over each member of theresulting subdivision of P .From [Loc16, lemma 5.2], for S1, S2 â S, S1 6= S2, and the points (x, y)such thatuS1(x, y) = uS2(x, y) = maxSâSuS(x, y)lie along a line as the functions uS1 and uS2 are linear. The quadratic andrational functions are the results of subproblems for h â Eâ²(P ), (a, b) âint(Dh) and the optimal solutions that could be shared by two subproblemswould always lie at the boundaries of the convex edge, that is, when (a, b) âDâh , D+h and over these regions the solution expressions are linear. Hence ifthe intersecting line passes through the polyhedral set over which both aredefined, it would lead to a further subdivision of that region again. However,the subdivision obtained after finding the maximum of all expressions wouldstill retain the polyhedral form.293.6. Algorithmic DesignSo finally, the closed convex envelope would again be characterized bya polyhedral subdivision such that over each member of the subdivision ithas one of the three functional forms discussed in Section 3.4.2.3.6 Algorithmic DesignOur algorithm is based on the same design as discussed in Section 3.4.1.An algorithm to compute the closed convex envelope using the above methodwould be quite complex in nature as it requires generating and potentiallysolving an exponential number of problems with sufficient floating pointaccuracy. Even a small floating point error could make an infeasible sub-problem feasible and vice versa.We also need to create a new module including appropriate data struc-tures and algorithms for all the 2D geometry manipulations required in theprocess. This involves methods for the map overlay problem which for Nlayers could be exponential in nature (see Section 3.3) and vertex enumer-ation problem to efficiently convert between halfspaces and vertices for apolytope.3.6.1 Input and Output Data StructuresThe data structure for the input quadratic function involves storing its6 coefficients in an array. For Q(x, y) = ax2 + bxy + cy2 + dx+ ey + f , westore the coefficients asQxy := [a b c d e f ].A linear function can be stored similarly by setting coefficients a = b = c = 0.The input polytope with N vertices is stored as an N Ã 2 matrix whereeach row contains a vertex and the rows are sorted in the clockwise direc-tion as shown in Figure 3.5. Since a 2D polytope is characterized by anedge-vertex (or vertex-edge) relationship, storing it requires storing all suchrelationships. One way would be to store these relationships as a graphwith adjacency lists. However, this will lead to redundant values in thedata structure. Another way would be to store the polytope as a set of lin-ear inequalities (halfspaces). This representation follows directly from thedefinition of polytope as being characterized by the intersection of a finitenumber of inequalities. These inequalities can be stored in a matrix, butin that case, obtaining the adjacent vertices and edges would require O (N)computations, whereas obtaining all the edge-vertex relationships would leadto O (N2) operations.303.6. Algorithmic Designa a x^2 + b xy + c y^2 + d x + ey + f22Q(x,y) = Qxy := a b c d e fA PolytopeP (Matrix)(x1,y1)(x2,y2)(x3,y3)(x4,y4)(x5,y5)x1 y1x yx2 y2x3 y3 x4 y4x5 y5Figure 3.5: Input Data StructureTo overcome this issue, we introduced the clockwise direction in ourdata structure. When vertices are stored in clockwise direction, obtainingthe neighborhood information is just a constant time operation as adjacentvertices would be at index ±1, and the edges between them can be obtainedby computing the equation of the line passing through the two sharing ver-tices, which is again O (1). We only stored the vertices as for a convexpolytope, by using the clockwise direction, we can obtain all the halfspacesin O (N) time and vice-versa.Example 3.6. Let Q(x, y) = 2x2+xyâ10y2+5xâ17yâ120 defined over apolytope formed by vertex set V = {(3, 10), (9,â5), (9, 10), (â8, 1), (7,â8)}.Then its input data structure representation would beQxy := [2, 1,â10, 5,â17,â120], andP :=ï£®ï£¯ï£¯ï£¯ï£¯ï£°3 109 109 â57 â8â8 1ï£¹ï£ºï£ºï£ºï£ºï£» .313.6. Algorithmic DesignOutput Polytope(x1,y1)(x2,y2)(x3,y3)(x4,y4)(x5,y5)(x6,y6)ax^2 + bxy + cy^2 + dx + ey + fux + vy + wFunctional form = 22Solxy(1,:) Solrgn{1,1}x1 y1x2 y2x6 y6x5 y5a1b1d1c1e1f1u1v1w1Solrgn{2,1}x2 y2x6 y622Solxy(2,:)a2b2d2c2e2f2u2v2w2x3 y3Solrgn{3,1}x6 y6x5 y522Solxy(3,:)a3 b3 d3c3 e3 f3 u3 v3 w3x3 y3 x4 y4Figure 3.6: Output Data StructureWe follow the same pattern for our output data structure. Since theclosed convex envelope is a piecewise defined rational function over a poly-hedrally divided domain, the coefficients of each function are stored in anarray, and each polytope in a vertex matrix. So forr(x, y) =ax2 + bxy + cy2 + dx+ ey + fux+ vy + w, the function is stored asRxy := [a b c d e f u v w],and a polytope with M vertices as a MÃ2 vertex matrix sorted in clockwisedirection. For a closed convex envelope characterized by N pieces, we storeall of the functions in a N Ã 9 matrix, and its corresponding polyhedralregions in a N Ã 1 cell with each element representing a polytope as shownin Figure 3.6.Example 3.7. Let Q(x, y) = xy defined over polytope P = {(â5, 5), (1, 3),323.6. Algorithmic Design(â1, 0), (â5,â4)}, then its closed convex envelope is given byconvQ,P (x, y) =ï£±ï£´ï£´ï£²ï£´ï£´ï£³5x2 â xy + 5y2 + 30xâ 30y + 25xâ y + 10 (x, y) â P115x2 â 3xy + 10y2 + 90xâ 65y + 753xâ 2y + 25 (x, y) â P2where P1 and P2 are triangles with vertices {(â5, 5), (â1, 0), (â5,â4)} and{(â5, 5), (1, 3), (â1, 0)}. So its output data structure representation isfxps(1,:) = [5,â1, 5, 30,â30, 25, 1,â1, 10],fxps(2,:) = [15,â3, 10, 90,â65, 75, 3,â2, 25],andrgns{1,1} =ï£®ï£°â5 5â1 0â5 â4ï£¹ï£» , rgns{2,1} =ï£®ï£°â5 51 3â1 0ï£¹ï£» .3.6.2 Main AlgorithmOur main algorithm is presented in Algorithm 1. It takes two inputarguments, Qxy and Pset, following the data structures discussed in Sec-tion 3.6.1. We start with verifying if the input quadratic function has anindefinite Hessian and then proceed onto transforming it into bilinear termusing the linear tranformations discussed in Section 3.1. We use the trans-formation matrix to obtain new transformed polytope points. After that,we separate the convex edges and the remaining vertices, and store themin the variable vertexSet. The vertexSet is a cell with two elements wherevertexSet{1, 1} contains the edges as rows in the form [v1, v2] with v1 and v2as the end points of the edge with v1x â¤ v2x , while vertexSet{2, 1} containsthe remaining vertices as rows.Example 3.8. Let Qxy = xy defined over polytope with Pset = [â5, 5; 1, 3;â1, 0;â5,â4]. Now the function is strictly convex over edges {(1, 3), (â1, 0)}and {(â1, 0), (â5,â4)}. SovertexSet{1, 1} =[â1 0 1 3â5 â4 â1 0], andvertexSet{2, 1} = [â5, 5].We use the vertexSet to compute the η functions according to convexedges and remaining vertices case, and store them in the variables convex-Eta and nonconvexEta respectively. Since η functions are either linear or333.6. Algorithmic Designquadratic (see Section 3.2.2) in a, b, we store the linear term xva + yvb + cas [xv yv c], and quadratic â(a+mjbâ qj)24mjâ bqj as [â mj qj ]. The âsignifies the quadratic form of η function.The element convexEta{1, 1} is a 3 Ã |Eâ²(P )| matrix containing all theη functions for the convex edges, while convexEta{2, 1} stores all the corre-sponding regions over which they are defined. Since the regions are eitherDâj = {(a, b) : a + bmj â¤ u1}, D+j = {(a, b) : l1 â¤ a + bmj â¤ +â} orDj = {(a, b) : l2 â¤ a + mj â¤ u2}, they are stored as [ââ, 1, mj , u1],[l1, 1, mj , +â] and [l2, 1, mj , u2] respectively.We follow the same pattern for nonconvexEta, but since they are definedover the full domain, their regions are stored as [ââ, 0, 0, +â].Example 3.9. For the function and polytope in example 3.8, the η functionsareη1(a, b) =ï£±ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£³a a+3b2â¤ â32â(a+3b2â 32)24â 3b2â32â¤ a+ 3b2â¤ 92âaâ 3b+ 3 92â¤ a+ 3b2,η2(a, b) =ï£±ï£´ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£´ï£³0 a+ b â¤ â9â(a+b2)22â9 â¤ a+ b â¤ â1â2aâ b+ 2 â1 â¤ a+ b,andη3(a, b) = 5aâ 5bâ 25.So they are stored asconvexEta{1, 1} =ï£®ï£¯ï£¯ï£¯ï£¯ï£¯ï£¯ï£°1 0 0â 1.5 1.5â1 â3 35 4 20â 1 11 0 0ï£¹ï£ºï£ºï£ºï£ºï£ºï£ºï£» ,343.6. Algorithmic DesignconvexEta{2, 1} =ï£®ï£¯ï£¯ï£¯ï£¯ï£¯ï£¯ï£°ââ 1 1.5 â1.5â1.5 1 1.5 4.54.5 1 1.5 âââ 1 1 â9â9 1 1 â1â1 1 1 âï£¹ï£ºï£ºï£ºï£ºï£ºï£ºï£» ,andnonconvexEta{1, 1} = [5,â5,â25],nonconvexEta{2, 1} = [ââ, 0, 0,â].The convexEta and nonconvexEta are then used to generate all the pos-sible subproblems. Since the subproblems are defined according to the sub-regions, first we solve the map overlay problem for all η functions and obtainall the subregions Srs. Then we use it to generate all the subproblems andstore them in the variable subprobs. Each subproblem, sp, is again storedas a cell with two elements. The sp{1, 1} contains the η functions, whilesp{2, 1} contains the subregions in the same form as convexEta.All the subproblems are solved using solve all subprobs() subroutinewhich follows the same approach as presented in Section 3.4.1. A detaileddescription of the underlying algorithm and data structures is covered inSection 3.6.3. The solutions of all the subproblems are then stored in thevariables pxps and pcdns where pxps is a cell containing all the solutionexpressions as cells while pcdns stores all the solution regions as a cell ofcells for the respective subproblems. Since the solution expressions are inthe rational form, they are stored as a 1 Ã 9 array in the cells of pxps. Onthe other hand, each cell of pcdns stores the regions as a set of inequalitiesin the form [yf,m, c, gf ] where yf and gf are flags for the y term and ââ¥âsign respectively. This form for storing halfspaces is discussed in Section 3.7.After that we have computed all the solutions, we solve the map overlayproblem to obtain the polyhedral subdivision of our transformed polytopenPset, and the maximum of all the expression over each member of subdivi-sion. We use pxps, pcdns and nPset for this purpose and store the respectiveexpressions and regions in nfxps and nrgns respectively. Then we use thelinear transformation matrix to revert back the transformations and obtainthe closed convex envelope with expressions stored in fxps and regions inrgns respectively. The data structure for fxps and rgns is covered in detailin Section 3.6.1.353.6. Algorithmic DesignAlgorithm 1 Computation of closed convex envelope for 2D Quadraticfunctions1: function cvxEnv2d(Qxy,Pset)2: qxp = Qxy(1:3);3: [isindefinite, isConvex] = hasIndefinite EigenVal(qxp);4: if isindefinite then5: quadxp = [0,1,0,0,0,0]; . Biliear form6: [Tm, lambda1 ] = transform quadxp to xy( Qxy );7: nPset = polytope points clockwise sort(Pset*Tm);8: vertexSet = find sconvex edgesn( quadxp, nPset);9: [convexEta, nonconvexEta] = compute eta( quadxp, vertexSet );10: if isempty(convexEta{1,1}) then11: Srs = compute subregions( convexEta{2,1} );12: ri = obtain subregion indices( Srs );13: subprobs = generate subproblems( convex-Eta,nonconvexEta,ri );14: else15: subprobs{1,1}{1,1} = nonconvexEta{1,1};16: subprobs{1,1}{2,1} = nonconvexEta{2,1}(1,:);17: end if18: [pxps, pcdns] = solve all subprobs( subprobs );19: if isempty(convexEta{1,1}) then20: [nfxps, nrgns, prgns] = compute polyhed subdiv( nPset,pcdns, pxps );21: else22: nfxps = pxps; nrgns = pcdns;23: end if[fxps, rgns] = revert transformations from xy( nfxps, nrgns, Tm,lambda1, Qxy, Pset);24: else25: if isConvex then26: fxps = [Qxy, 0, 0, 1]; rgns{1,1} = Pset; . Convex case27: else28: print âFunction is Concave.â . Function is concave!29: end if30: end if31: return fxps, rgns32: end function363.6. Algorithmic Design3.6.3 Solving the subproblemsThe algorithm to solve subproblems over a subregion, at its core, followsthe same approach discussed in section 3.4.1. In this section, we discuss thedata structures and design decisions made for solving the subproblems inaccordance with Algorithm 1.In Algorithm 2, each subproblem sp is a cell of two elements, both ofwhich are matrices. The matrix sp{1,1} stores all the coefficients of the ηfunctions while sp{2,1} contains their respective subregions Sr. For eachη pair, we solve the equality in cvxHull2d equality solver() and store therelation in the variable relab. The expressions for a and b are stored ina cell with relab{1,1} and relab{2,1} containing their respective relations.The element relab{1,1} is a vector(or matrix) of size 3(or 2 Ã 3) containingthe coefficients of the quadratic or linear relation(s) in Z or b. The elementrelab{2,1} is empty if a is defined in terms of b, otherwise it is a vector ofsize 3 containing the quadratic relation in Z, similar to relab{1,1}.Using the relation, we reduce the subregion Sr to one-dimension. Ifthe one-dimension region is feasible( nonempty ), we solve the inequalityconstraint for the remaining ηj . When all the inequalities are feasible, weformulate the one-dimension optimization problem in Z or b knowing that wewill have a valid solution too. The problem is computed using compute OP()and is stored in the variable mOP as a matrix of size 4Ã3. Each row containsthe coefficients of x,y and the constant term which are the coefficients of thepowers of the optimization variable u respectively(see Section 3.4.3). Row 1contains the coefficients of quadratic, row 2 of linear, and row 3 the constantterm in the optimization variable. Row 4 contains the case number of theproblem for our solver function OP Solver().Example 3.10. For an optimization problem of formmax(x36â y36+518)Z2 + 2(x4+y4)Z +11x4â 11y4+52â 9 â¤ Z â¤ â1defined in Z would be stored asmOP =ï£®ï£¯ï£¯ï£°1/36 â1/36 5/181/4 1/4 011/4 â11/4 5/20 1 0ï£¹ï£ºï£ºï£» .The fourth row signifies that the optimization problem comes from thequadratic-linear case, and is defined in Z.373.6. Algorithmic DesignAlgorithm 2 Algorithm for solving subproblems over one subregion1: function solve subprobs over subregion(sp)2: ae = sp{1,1}; ac = sp{2,1}; sc = 1;3: for i=1 to size(ae) do4: for j=i+1 to size(ae) do5: eh = ae(i,:); ew=ae(j,:);6: [ relab, isfeasibleRelab] = cvxHull2d equality solver(eh,ew);7: if isfeasibleRelab then8: for ii=1:size(relab) do9: ab = relab{ii,1};10: [Sr1d, isfeasibleSr1d] = compute Sr 1d(ac,ab);11: if isfeasibleSr1d then12: [feasible region, isfeasible] = cvx-Hull2d inequality solver(eh, ew, Sr1d, ab, ae);13: if isfeasible then14: [mOP] = compute OP( eh, ew, ab );15: [ txps, tcdns ] = OP Solver( mOP, feasi-ble region );16: solxp{sc,1} = txps; solcdn{sc,1} = tcdns;17: sc = sc + 1;18: end if19: else20: print âThe subproblem for eh,ew is infeasible.â21: end if22: end for23: end if24: end for25: end for26: end function383.6. Algorithmic DesignSimilarly an example of optimization problem resulting from solving alinear-linear equality in b would bemax(5x4+ y â 54)b+(25x4â 254)â 299â¤ b â¤ â3111,and it would be stored asmOP =ï£®ï£¯ï£¯ï£°0 0 05/4 1 5/425/4 0 25/40 1 1ï£¹ï£ºï£ºï£»The feasible regions are stored in the feasible region variable.We use the solver function OP Solver() to solve the optimization prob-lem mOP and obtain the respective solution expressions and regions. Theexpressions are stored in txps as a cell of 1Ã9 vectors and the regions in tcdnsas a cell of matrices containing the region in the form l â¤ ax+ by+ c â¤ u as[l, a, b, c, u, ef ] where the ef flag differentiates between ââ¤â and â<â with 1for the former and 0 for the latter case.Example 3.11. For the Optimization problem defined in Z in example 3.10,the solution would be in two regions, that is, whenx36â y36+518> 0, thefunction would beï£±ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£³5x2 â xy + 30x+ 5y2 â 30y + 25xâ y + 10 xâ y + 10 < x+ y < â19(xâ y + 10)â4xâ 5y â 20 xâ y + 10 â¥ x+ y20xâ 29y + 209â19(xâ y + 10) â¤ x+ y,and whenx36â y36+518â¤ 0, the function would bemax{â4xâ 5y â 20, 20xâ 29y + 209}.393.7. Map Overlay problemSo they would be stored in the variables txps and tcdns astxps{1,1} = [5/36,â1/36, 5/36, 5/6,â5/6, 25/36, 1/36,â1/36, 5/18],tcdns{1,1} =ï£®ï£° 0 1/36 â1/36 5/18 â 00 1/2 0 5/2 â 0ââ 5/18 2/9 5/18 0 0ï£¹ï£» ,txps{2,1} = [0, 0, 0,â4,â5,â20, 0, 0, 1],tcdns{2,1} =[0 1/36 â1/36 5/18 â 0ââ 1/2 0 5/2 0 1],txps{3,1} = [0, 0, 0, 20/9,â29/9, 20/9, 0, 0, 1],tcdns{3,1} =[0 1/36 â1/36 5/18 â 00 5/18 2/9 5/18 â 1],txps{4,1} =[0, 0, 0,â4,â5,â20, 0, 0, 10, 0, 0, 20/9,â29/9, 20/9, 0, 0, 1],tcdns{4,1} = [ââ, 1/36,â1/36, 5/18, 0, 1].More than one expression in an expression matrix signifies the maximum ofall, for example in txps{4,1} the solution expression is defined in two linearexpressions. Also, when we have more than one expression, we may havefurther subdivision of the region as well.3.7 Map Overlay problemMap overlay is one of the most commonly occurring problem when work-ing with piecewise functions in more than one-dimension. As simple asadding two piecewise functions requires solving a map overlay problem ontheir domains. The data structures used to store the domains play a sig-nificant role in simplifying the problem a bit. The map overlay problem isexponential in nature, that is, for N layers with m subdivisions each, theresulting subdivision can have mN subregions. The type of domain divisionof the functions can further add some additional complexity to the problem.For example, to solve the map overlay problem for two parabolic subdivi-sions, not only do we need to enumerate all the possible intersections butalso add some more divisions to the solution to properly store some sub-regions. However, to obtain the closed convex envelope of a 2D quadraticfunction, we only need to solve the problem for polyhedral subdivisions.Next, we design an efficient data structure to reduce the complexity of theproblem in our module.403.7. Map Overlay problemy â¥ mx + cy â¤ mx + cy =mx + cFigure 3.7: Space division by line y = mx+ c3.7.1 Data structures for polyhedral regionThe most common form for representing an inequality is ax + by â¤ cwhere a, b, c â R. The form can be stored in an array as [a b c]. For u âRn, it is cTu â¤ d where c â Rn and d â R, and the data structure canbe expanded to size n + 1 to store inequalities of this form. This datastructure performs well as it is intuitive too, however the representation ofany unique inequality, in this case, is not unique. For example, ax+ by â¤ cand N(ax+ by â¤ c) âN â R>0 are two different representations of the sameinequality. To force unicity of the form, it can be normalized by dividingthe inequality with the coefficient of either x or y term. But in that case,not all inequalities can be represented by this form, we would need ââ¥â signtoo.In our algorithm, we came across the map overlay problem twice. First,when generating all the possible subproblems in a â b space and second,when computing the polyhedral subdivision for the closed convex envelopein the given polytope. For generating the subproblems, we solve the mapoverlay problem for η functions with 3 polyhedral regions separated by twoparallel hyperplanes. We stored the subregions l â¤ a+bmj â¤ u as [l, 1,mj , u]while regions a + mj â¤ u and l â¤ a + mjb are stored as [ââ, 1,mj , u] and[l, 1,mj ,â] respectively. As this is a normalized form too, there always willbe a unique representation for each inequality, however it cannot be used tostore any general inequality in two variables.413.7. Map Overlay problemTo do so, we formulate it as either y â¤ mx+ c or y â¥ mx+ c and thenstore it as [yf,m, c, gf ]. The motivation for this representation comes fromthe general equation of a line y = mx + c. This hyperplane divides the R2space into two halfspaces as shown in Figure 3.7. This form is normalizedas well because by forcing the coefficient of y to 1, we make sure that therecan only be one representation of an inequality. Any other representationobtained by multiplying a number N â R>0 will still have the same form.When the coefficient of y is zero, we make the coefficient of x as 1, therebymaintaining the unique representation.Now to store this form, we need four values [yf,m, c, gf ] compared toonly three in [a, b, c] because for uniqueness we need a flag for the ââ¥âsign. The variable yf is a flag for the y term, if it is 1 then the inequalityis defined in terms of x and y, otherwise it would only be in x. The termsm and c are corresponding terms from the inequality with m forced to be 1whenever yf = 0. The gf flag when 1 represents ââ¥â and when 0 representsthe ââ¤â case.Example 3.12. For y â¥ 2x + 3, y â¤ 3x + 4, 0 â¥ x + 2, y â¥ 1, their datastructure arrays would be [1, 2, 3, 1], [1, 3, 4, 0], [0, 1, 2, 1], and [1, 0, 1, 1] .Using the four terms in our data structure, we can uniquely representany general inequality with the least amount of space possible (two of themare booleans). The other advantage is that it is easier to substitute thevalue of y in our algorithms. In our module, we also have a subroutineconvert inequalities() which can convert inequalities of form l â¤ ax+by+c â¤u to [yf,m, c, gf ] form.3.7.2 Sorting vertices in clockwise directionAs discussed before in section 3.6.1, storing points in clockwise direc-tion significantly reduces the amount of computation required to obtain theneighborhood information for a vertex or edge. Apart from being an inputdata structure, our algorithms heavily rely on the computation of the newpolytopes, so we need a method to sort the vertices in clockwise directionas well.The algorithm to sort polytope vertices in clockwise direction is pre-sented in Algorithm 3. The idea has been taken from the working of a clockitself. For a given set of vertices, we first compute the centroid. Since thecentroid always would lie inside the polytope (as it is convex), we use thispoint as the center of the clock.423.7. Map Overlay problemAlgorithm 3 Sorting vertices of a polytope in clockwise direction1: function polytope points clockwise sort(P)2: cen = mean(P); ps = []; ns = [];3: for i=1 to size(P) do4: ts = compute slope(P(i,:), cen);5: if ts â¥ 0 or isinf(ts) then6: ps = [ps; [P(i,:), abs(ts)]]; . positive slopes or Â±â7: else8: ns = [ns; [P(i,:), ts]]; . negative slopes9: end if10: end for11: [ q1p, q3p ] = polytope points sort Q13(ps, cen); . Sorts accordingto slopes in quadrant 1 & 312: [ q2p, q4p ] = polytope points sort Q24(ns, cen); . Sorts accordingto slopes in quadrant 2 & 413: cwp = [ q2p; q1p; q4p; q3p ];14: return cwp15: end functionWe divide the R2 space into four quadrants using centroid as the origin.Then we compute the slope of each line connecting the centroid with thevertices. Now the slopes for points in quadrants 1 and 3 will be positivewhile in quadrants 2 and 4 will be negative. The special cases of slopes as0 or Â±â have been accommodated in the quadrants 1 and 3 case itself. Wesort the points in quadrants 1 and 3, and quadrants 2 and 4 according totheir slopes in the clockwise direction (opposite of how we measure anglesconventionally) and obtain the sorted points in all the 4 quadrants. Theclockwise sorted points are stored in the variable cwp in the order: quadrant2, 1, 4 and 3 respectively. It does not matter how the ordering begins aslong as they follow the clockwise order.Example 3.13. For a polytope P = {(8, 10), (10, 7), (7, 8), (1, 8), (9, 4)}, westart with computing the centroid cen = (7, 7.4). Next, we compute theslope of all the lines connecting the points with the centroid as433.7. Map Overlay problemcen(2)â P (i, 2)cen(1)â P (i, 1) to getsv =ï£®ï£¯ï£¯ï£¯ï£¯ï£°13/5â2/15âââ1/10â17/10ï£¹ï£ºï£ºï£ºï£ºï£» .Assuming cen as the origin, we divide the polytope points according topositive slopes (including Â±â) and negatives slopes, that is,psp =[8 10 13/57 8 ââ], nsp =ï£®ï£°10 7 â2/151 8 â1/109 4 â17/10ï£¹ï£» .In the algorithm, we convert ââ to â, so that psp contains all the pointswith positive slopes. This conversion has been properly accommodated inpolytope points Q1 3(). After that, we divide the points with positive slopesin quadrants 1 and 3, and with negative slopes in quadrants 2 and 4. So thepoints and slopes combination becomesq1 =[8 10 13/57 8 ââ], q3 =[ ]andq2 =[1 8 â1/10], q4 = [10 7 â2/159 4 â17/10].Finally, after sorting the points in each quadrant according to arctan of theirslopes, in decreasing order, we getq1p =[7 88 10], q3p =[ ]andq2p =[1 8], q4p =[10 79 4].So the polytope points in the clockwise order (according to algorithm default443.7. Map Overlay problemorder is [q2p; q1p; q4p; q3p]) arecwp =ï£®ï£¯ï£¯ï£¯ï£¯ï£°1 87 88 1010 79 4ï£¹ï£ºï£ºï£ºï£ºï£» .The time complexity of this algorithm is bounded by the time requiredto sort the points by decreasing slopes in each quadrant. So for N vertices,the worst case time complexity would be O (N logN). Although it takesO (N logN) time to sort the points, it would still perform better compared tounsorted points as it would take O (1) operation to obtain the neighborhoodinformation which would have taken O (N) operations in the latter case.Also, the transformation between vertices and halfspaces representationswould only take O (N) computations, which otherwise would have beenO (N2) with unsorted points.3.7.3 Map Overlay AlgorithmSince for a map overlay problem for N layers with m subdivisions each,the resulting region could have mN number of subdivisions in the worst case,any algorithm aiming at finding the solution cannot perform better than abrute force method that enumerates all the possible subregions. However,still a constructive algorithm can be designed which builds the solution byadding regions on top of another and gets rid of the infeasible subregions ateach step. An example with three layers have been shown in Figures 3.1-3.4.The main algorithm for map overlay problem is presented in Algorithm 4.Its implementation in the module has API compute polyhed subdiv() as alongwith computing the subregions, we need to store the corresponding functionexpressions as well. However, the subroutine map overlay prob() presentedhere contains the main concept of the additive computation of subregions.At each iteration, the computation of the solution subregions is per-formed in compute intersection subregions(). For two layers Sr1 and Sr2,we first enumerate all the possible subregions pair, compute the intersectionbetween each, and then remove the infeasible ones. If Sr1 and Sr2 have n1and n2 subregions, then each call to this method takes O (n1 Ã n2) time.Now, at the subregion level, for two polyhedral sets sr1 and sr2, theresulting subregion is the intersection of their inequality sets. In our case,all the subregions are bounded polyhedral sets (as they are a member of thesubdivision of the input polytope). So we enumerate the inequality pairs453.7. Map Overlay problemAlgorithm 4 Constructive algorithm for map overlay problem for polytopes1: function map overlay prob(rgns)2: Srs = {}; . Initialize Srs as empty cell3: N = size(rgns); . Number of regions4: if N==1 then . If only one region5: Srs = rgns;6: else7: for i=1 to N-1 do8: if i==1 then9: Srs=compute intersection subregions(rgns{1,1},rgns{2,1});10: else11: Srs = compute intersection subregions(Srs, rgns{i+1,1});12: end if13: end for14: end if15: return Srs16: end function17:18: function compute intersection subregions(Sr1,Sr2)19: nSrs = ;20: cnt = 1; . Counter for new subregions21: for i=1 to size(Sr1) do22: for j=1 to size(Sr2) do23: tmpsr=compute intersection of two subregions(Sr1{i,1},Sr2{j,1});24: tmpPoly = compute polytope from regions(tmpsr);25: if size(tmpPoly) > 2 then26: nSrscnt,1 = compute polytope region(tmpPoly);27: cnt = cnt + 1;28: end if29: end for30: end for31: return nSrs32: end function463.7. Map Overlay problemAlgorithm 5 Computing intersection of two subregions1: function compute intersection of two subregions(sr1,sr2)2: isIntersectionfeasible = true;3: for i=1 to size(sr1) do4: for j=1 to size(sr2) do5: isfeasible = intersection of two halfspaces( sr1(i,:),sr2(j,:));6: if !isfeasible then7: isIntersectionfeasible = false;8: break;9: end if10: end for11: if !isfeasible then12: break;13: end if14: end for15: if isIntersectionfeasible then16: nsr = reduce to required inequalities([sr1;sr2]);17: if isempty(nsr) then18: isIntersectionfeasible = false;19: end if20: else21: nsr = [];22: end if23: return nsr24: end function473.7. Map Overlay problemFigure 3.8: Vertex Enumeration Problemfrom both sets, then verify if they are feasible, if at least one pair is not, weset the isfeasible flag to false and return an empty set. When all the pairsare feasible, we reduce the number of inequalities to the only active onesand return it using nsr.The problem of reducing the inequalities to only the active inequalities,a variant of the Vertex Enumeration Problem, is known to be an NP-hardproblem for unbounded polyhedral sets [KBB+09, Lin86], however, we onlyhave to deal with polytopes, especially 2D polytopes, as the subregions area member of the polyhedral subdivision of the input polytope only.The method of computing the vertices is the same as the intersectionof two subregions. For each inequality pair, we compute their intersectionpoint and verify if it satisfies all other inequalities. Since for each pair,we are verifying the point with every other inequality, the computation ofpolytope takes O (N3) operations for N number of inequalities.The closed convex envelope of a 2D nonconvex quadratic function is apiecewise rational function defined over a polyhedral subdivision of the inputpolytope. After computing the convex envelope, we proceed to compute theconjugate of each piece of the closed convex envelope in the next chapter.48Chapter 4Conjugate of a class ofConvex Bivariate RationalFunctions over a polytopeThe closed convex envelope of a bivariate nonconvex quadratic functiondefined over a polytope is a piecewise rational function of form (3.18) witha polyhedral subdivision. The conjugate of the convex piecewise rationalfunction, say R, is then computed asRâ(s) = maxjâ{1,...,N}Râj (s)where N is the total number of pieces and Râj (s) is the conjugate of the jthpiece.In this chapter, we discuss the computation of the conjugate of a convexrational function defined over a polytope. We start with computing the sub-differentials corresponding to each entity of the polytope and then computethe expressions over them to finally obtain the conjugate. All the functionsinvolved here are in convRn.Lemma 4.1. For x â Rn, a convex set P â Rn, and the indicator functionIP (x) ={0 x â P+â otherwise,the subdifferential set âIP (x) = NP (x), where NP (x) is the normal cone ofP at point x.Proof. The subdifferential of IP (x) isâIP (x) = {s â Rn : IP (y) â¥ IP (x) + ãs, y â xã, ây â P}= {s â Rn : 0 â¥ 0 + ãs, y â xã, ây â P}= {s â Rn : ãs, y â xã â¤ 0, ây â P}= NP (x).494.1. Subdifferentials in the interior of the polytopeLemma 4.2. For a function g : Rn â R âª {+â} with dom(g) = Rn, apolyhedral set P â Rn, let f(x) = g(x)+IP (x), then âf(x) = âg(x)+NP (x)andâf(x) =ï£±ï£´ï£²ï£´ï£³âg(x) x â int(P )âg(x) +NP (x) x â bndry(P )Ï otherwise.Proof. For f(x) = g(x) + IP (x), the subdifferential set âf(x) = â(g(x) +IP (x)). Since dom(g) = Rn, so dom(g) â© P = P . Hence from [HUL13, p.114, Corollary 3.1.2], for all x â Pâf(x) = â(g(x) + IP (x))= âg(x) + âIP (x)= âg(x) +NP (x) [using Lemma 4.1].Now for any x â int(P ), NP (x) = {0}. So âf(x) = âg(x).In our case, the function g will either be a convex rational, quadraticor linear function of forms (3.18)-(3.20) respectively. In addition, g is theclosed convex envelope of a bivariate nonconvex function Q satisfying As-sumption 3.2.1.4.1 Subdifferentials in the interior of thepolytopeProposition 4.3. For a bivariate convex quadratic nonlinear function q ofform (3.19), there exists α10, α01, α00 â R such thatâxâR2 âq(x) = {s :Cq(s) = 0} where Cq(s) = α10s1 + α01s2 + α00.Proof. From form (3.19), let q(x) =ξ21(x)ξ20+ ξ0(x) be a convex quadraticfunction with ξ1(x) = ξ11x1 + ξ12x2 + ξ10, ξ0(x) = ξ01x1 + ξ02x2 + ξ00, ξ20 >0, and ξ11 or ξ10 6= 0. The parameter ξ20 is included to match αij withFormula (3.19).Since q â C1, for any x â R2âq(x) = {s : s = âq(x)} .504.1. Subdifferentials in the interior of the polytopeThe gradient s = âq(x) is then computed ass1 =2ξ11ξ20t+ ξ01s2 =2ξ12ξ20t+ ξ02,(4.1)where t = ξ11x1 + ξ12x2 + ξ10.The Equation (4.1) represents the parametric equation of a line in theplane and by eliminating t, we getâÎ¾12s1 + ξ11s2 + ξ01ξ12 â ξ02ξ11 = 0which is the Cartesian equation of our line. Let Cq(s) = α10s1 +α01s2 +α00with α10 = âÎ¾12, α01 = ξ11, and α00 = ξ01ξ12 â ξ02ξ11. Then for allx â R2, âq(x) is contained in the line Cq(s) = 0, that is,âxâR2âq(x) â {s : Cq(s) = 0}. (4.2)Now consider a point sq that satisfies Cq(sq) = 0, so it must satisfy theparametric equation of the same line as well. Hence we also have{s : Cq(s) = 0} ââxâR2âq(x). (4.3)So by using (4.2) and (4.3) we getâxâR2âq(x) = {s : Cq(s) = 0}.Example 4.4. For a convex quadratic function q(x) =14(x1+x2)2+18(x1âx2), the gradient s = âq(x) iss1 =t2+18, s2 =t2â 18(4.4)where t = x1 + x2. Upon solving the above equations, we getâ4s1 + 4s2 + 1 = 0.514.1. Subdifferentials in the interior of the polytopePrimal domain (x1, x2) Dual domain (s1, s2)x1 + x2 = 0x1 + x2 â 1 = 0x1 + x2 + 1 = 0â4s 1+4s 2+ 1= 0(18 , 18)(58 , 38)(â 38 , â 58)Figure 4.1: Primal and dual mapping for quadratic function in Example 4.4over lines parallel to x1 + x2 = 0.For any c â R, let t â c = 0, that is, x1 + x2 â c = 0 be a line parallel tox1 + x2 = 0. When c = 0, both lines coincide. Now (4.4) can be written ass1 =c2+18, s2 =c2â 18.So the subdifferentials along each dimension are a linear function of theintercepts of the lines parallel to x1 + x2 = 0. This is well illustrated inFigure 4.1. The whole primal domain (shown in grey) maps to the lineâ4s1 + 4s2 + 1 = 0 (also grey) in the dual domain. The line to pointmapping is shown by the three lines x1 + x2 = 0, x1 + x2 â 1 = 0 andx1 +x2 + 1 = 0 of color red, green and blue respectively, to the points of thesame color on the line â4s1 + 4s2 + 1 = 0 in the dual domain.Corollary 4.5. For a linear function l of form (3.20), and any x â R2, thesubdifferential set âl(x) is a singleton.Proof. By setting ξ11 = ξ12 = 0 in Proposition 4.3, the Equation (4.1)reduces tos1 = ξ01 and s2 = ξ02.Hence âl(x) = {(ξ01, ξ02)} is a singleton set.524.1. Subdifferentials in the interior of the polytopeCorollary 4.6. For a bivariate convex quadratic function q of form (3.19),and a polytope P â R2, let f(x) = q(x) + IP (x), then for all x â int(P ), thesetâxâint(P ) âf(x) is contained in a line segment.Proof. Since P â R2 is a compact set and q â C1, âq(P ) is compact aswell. Using Proposition 4.3,âxâP âq(x) â L where L â R2 is a line. Sinces is defined by the parametric equation of a line and P is convex (and sois connected),âxâP âq(x) = âq(P ) is a segment andâxâint(P ) âf(x) ââxâP âq(x) is contained in the same line segment as well.Proposition 4.7. For a bivariate rational function r of form (3.18), thereexists αij such that the setâxâdom(r) âr(x) = {s : Cr(s) = 0}, where Cr(s) =α11s21 + α12s1s2 + α22s22 + α10s1 + α02s2 + α00 and Cr(s) = 0, is a paraboliccurve (as per Definition 2.20).Proof. From form (3.18), let r(x) =ξ21(x)ξ2(x)+ ξ0(x) be a bivariate rationalfunction with ξ1(x) = ξ11x1 + ξ12x2 + ξ10, ξ2(x) = ξ21x1 + ξ22x2 + ξ20and ξ0(x) = ξ01x1 + ξ02x2 + ξ00. Since r is differentiable everywhere indom(r) = R2/{z : ξ21z1 + ξ22z2 + ξ20 = 0}, for any x â dom(r)âr(x) = {s : s = âr(x)}.Next we compute the gradient s = âr(x) ass1 = 2Î¾11tâ ξ21t2 + ξ01s2 = 2Î¾12tâ ξ22t2 + ξ02(4.5)where t =ξ11x1 + ξ12x2 + ξ10ξ21x1 + ξ22x2 + ξ20. Equation (4.5) represents the parametricequation of a conic section and by eliminating t, we getα11s21 + α12s1s2 + α22s22 + α10s1 + α02s2 + α00 = 0 (4.6)where α11 = ξ221ξ222, α12 = â2Î¾321Î¾22, α22 = ξ421 and other αij are functionsof the coefficients of r.Now sinceα212 â 4α11α22 = (â2Î¾321Î¾22)2 â 4ξ621ξ222 = 0,the conic section represented by (4.6) is a parabola.Let Cr(s) = α11s21 + α12s1s2 + α22s22 + α10s1 + α02s2 + α00, then for allx â dom(r), âr(x) is contained in the parabolic curve Cr(s) = 0, so534.1. Subdifferentials in the interior of the polytopeâxâdom(r)âr(x) â {s : Cr(s) = 0}. (4.7)Now any point sr that satisfies Cr(sr) = 0, satisfies the parametricequation (4.5) as well. So we also have{s : Cr(s) = 0} ââdom(r)âr(x). (4.8)Hence by using (4.7) and (4.8), we getâxâdom(r)âr(x) = {s : Cr(s) = 0}.Example 4.8. For a rational function r(x) =x22x2 â x1 + 1 and any x âdom(r), the gradient s = âr(x) iss1 = t2, s2 = 2tâ t2where t =x2x2 â x1 + 1. After eliminating t, the equation of the parabolacomes out to bes21 + 2s1s2 + s22 â 4s1 = 0.So from Proposition 4.7,âxâdom(r)âr(x) = {s : s21 + 2s1s2 + s22 â 4s1 = 0}.Unlike the quadratic function in Example 4.4, the mapping between primaland dual domain is difficult to visualize for rational functions.Corollary 4.9. For a bivariate rational function r, of form (3.18), a poly-tope P , let f(x) = r(x) + IP (x), then for all x â int(P ), the setâxâint(P )âf(x) is contained inside a parabolic arc.Proof. Using the same arguments from Corollary 4.6,âxâint(P ) âf(x) ââxâint(P ) âr(x) ââxâdom(r){âr(x)}. From Proposition 4.7,âxâint(P ) âr(x)â P where P â R2 is a parabolic curve. Since s is defined by the para-metric equation of a parabola and P is connected, soâxâint(P ) âr(x) =âr(int(P )) â P is contained in a parabolic arc.544.2. Subdifferentials at the Vertices4.2 Subdifferentials at the VerticesLemma 4.10. For a function g â C1, a polytope P , and vertex v, let f(x) =g(x) + IP (x), then the subdifferential set âf(v) is an unbounded polyhedralset.Proof. Let v be a vertex of P , and vâ and v+ be the adjacent verticeswhen traversed P in the clockwise direction, then for all y â P , yl â Eâand yr â E+ where Eâ and E+ are the respective adjacent edges to v, thenormal cone at v isNP (v) = {s â R2 : ãs, y â vã â¤ 0, ây â P}= {sl â R2 : ãsl, yl â vã â¤ 0} â© {sr â R2 : ãsr, yr â vã â¤ 0}= {s â R2 : ãs, yl â vã â¤ 0, ãs, yr â vã â¤ 0},which is an unbounded polyhedral set. Now from Proposition 4.2, we haveâf(v) = âg(v) +NP (v)= âg(v) +NP (v)= {s+âg(v) â R2 : ãs, yl â vã â¤ 0, ãs, yr â vã â¤ 0}which is again an unbounded polyhedral set. It is illustrated as the arrowsin the intersection of Red and Green regions in Figure 4.2.Conjecture 4.11. For a bivariate rational function r of form (3.18), letf(x) = r(x) + IP (x) and v â P be a vertex of P such that ξ1(v) = ξ2(v) =0, then the subdifferential set âf(v) is a parabolic region (as per Defini-tion 2.21)(Note that r is not differentiable at v so Lemma 4.10 does notapply).4.3 Subdifferentials on the edgesLemma 4.12. For a bivariate function g â C1, a polytope P , and an edgeE = {x : x2 = mx1 + c, xl1 â¤ x1 â¤ xu1} between vertices xl and xu, letf(x) = g(x)+IP (x), thenâxâri(E) âf(x) =âxâri(E){s+âg(x) : s1+ms2 =0, s2 â¥ 0}.Proof. From Lemma 4.2, for all x â ri(E), âf(x) = âg(x) + NP (x). LetL(x) = x2 âmx1 â c be the expression of the line joining xl and xu suchthat P â {x : L(x) â¤ 0}. The case P â {x : L(x) â¥ 0} is analogous.554.3. Subdifferentials on the edgesFigure 4.2: Normal cone at a vertex v â P shown as the arrows in theintersection of Red and Green regionSince P â R2 is a polytope, for all x â ri(E), NP (x) = {s : s =Î»âL(x), λ â¥ 0} is the normal cone of P at x. SoNP (x) = {s : s = Î»âL(x), λ â¥ 0}= {s : s = Î»(âm, 1), λ â¥ 0}= {s : s1 = âmÎ», s2 = λ, λ â¥ 0}= {s : s1 = âms2, s2 â¥ 0} [by eliminating λ]= {s : s1 +ms2 = 0, s2 â¥ 0}.In the special case when E = {x : x1 = d, xl1 â¤ x1 â¤ xu1}, L(x) = x1 â d,soNP (x) = {s : s = λ(1, 0), λ â¥ 0}= {s : s2 = 0, s1 â¥ 0}.Now for any x â ri(E),âf(x) = âg(x) +NP (x)= âg(x) + {s : s1 +ms2 = 0, s2 â¥ 0}= {s+âg(x) : s1 +ms2 = 0, s2 â¥ 0}.(4.9)564.3. Subdifferentials on the edgesSo âxâri(E)âf(x) =âxâri(E){s+âg(x) : s1 +ms2 = 0, s2 â¥ 0}.Proposition 4.13. For a true bivariate convex quadratic function q ofform (3.19), a polytope P and an edge E = {x : x2 = mx1 + c, xl1 â¤x1 â¤ xu1} between vertices xl and xu, let f(x) = q(x) + IP (x), then the setâxâri(E) âf(x) is either an unbounded polyhedral set with nonempty interioror a ray.Proof. From Corollary 4.6, there exists some l, u such thatâxâri(E) âq(x) =âxâri(E){s : Cq(s) = 0, l1 â¤ s1 â¤ u1}. Now computation ofâxâri(E) âf(x)can be divided in the following two cases:Case 1 (l = u) From Corollary 4.6,âxâri(E)âq(x) =âxâri(E){s : Cq(s) = 0, l1 â¤ s1 â¤ u1}= {s : Cq(s) = 0, s = l = u = âq(x)} [since l = u]= uSo from Lemma 4.12,âxâri(E)âf(x) =âxâri(E){s+âq(x) : s1 +ms2 = 0, s2 â¥ 0}= {s+ u : s1 +ms2 = 0, s2 â¥ 0}= {s : s1 +ms2 â (u1 +mu2) = 0, s2 â¥ u2}(4.10)is a ray.Case 2 (l 6= u) From Corallary 4.6, âxâri(E) âq(x) = {s : Cq(s) = 0, l1 â¤s1 â¤ u1}. Now for any x â ri(E), âf(x) = {s : s1 + ms2 = 0, s2 â¥ 0}. Sowhen âq(x) = l, âf(x) = {s : s1 +ms2 â (l1 +ml2) = 0, s2 â¥ l2} and whenâq(x) = u, âf(x) = {s : s1 + ms2 â (u1 + mu2) = 0, s2 â¥ u2}. Now letâf(x) â {s : Cq(s) â¤ 0}, the case âf(x) â {s : Cq(s) â¥ 0} is analogous.Thenâxâri(E)âf(x) =âxâri(E){s+âq(x) : s1 +ms2 = 0, s2 â¥ 0}= {s : l1 +ml2 â¤ s1 +ms2 â¤ u1 +mu2, Cq(s) â¤ 0}(4.11)is an unbounded polyhedral set with nonempty interior.574.3. Subdifferentials on the edgesExample 4.14. For the quadratic function in Example 4.4, and polytopeP with vertices v1 =(12, 1), v2 =(34,34), v3 =(14,14)and v4 =(0,12),let f(x) = q(x) + IP (x).For an edge E1 ={x : x2 = âx1 + 32,12â¤ x1 â¤ 34}between vertices v1and v2, ââf(v1) = ââf(v2) =(78,58), as the edge is parallel to x1+x2 = 0.So in this case,âxâri(E) âf(x) = {s : s1 â s2 â 1/4 = 0, s2 â¥ 5/8} is a ray.For the edge E2 ={x : x2 = x1,14â¤ x1 â¤ 34}between vertices v2 andv3, âf(v2) =(78,58)and âf(v3) =(38,18). Also, from Example 4.4,Cq(s) = â4s1 + 4s2 + 1. Now âf(v2) = {s : s1 â s2 â 1/4 = 0, s2 â¥ 5/8}and âf(v2) = {s : s1 â s2 â 3/2 = 0, s2 â¥ 1/8}. Clearly, for any x â ri(E),âf(x) â {s : Cq(s) â¤ 0}. Soâxâri(E)âf(x) ={s :12â¤ s1 â s2 â¤ 32,â4s1 + 4s2 + 1 â¤ 0},which is an unbounded polyhedral set.Corollary 4.15. For a bivariate linear function l, a polytope P and anedge E = {x : x2 = mx1 + c, xl1 â¤ x1 â¤ xu1} between vertices xl and xu, letf(x) = l(x) + IP (x), then the subdifferential setâxâri(E) âf(x) is a ray.Proof. From Corollary 4.5, for all x â R2, âl(x) is a singleton {âl(x)}. Letc = âl(x). Then from Lemma 4.12,âxâri(E)âf(x) =âxâri(E){s+âl(x) : s1 +ms2 = 0, s2 â¥ 0}= {s : s1 +ms2 â (c1 +mc2) = 0, s2 â¥ c2}is a ray.Proposition 4.16. For a bivariate rational function r of form (3.18), apolytope P and an edge E = {x : x2 = mx1 + c, vâ1 â¤ x1 â¤ v+1 } betweenvertices vâ and v+, let f(x) = r(x) + IP (x), then the subdifferential setâxâri(E) âf(x) is either a parabolic region(see Definition 2.21) or a ray.Proof. From Corollary 4.9, there exists l, u â R2 such that âxâri(E) âr(x) =âxâri(E){s : Cr(s) = 0, l1 â¤ s1 â¤ u1}. So computingâxâri(E) âf(x) leads tothe following two cases:584.4. Structure of the conjugate domainCase 1 (l = u) Same as Proposition 4.13, Case 1.Case 2 (l 6= u) By setting g = r in Lemma 4.12, for any x â ri(E),âf(x) = {s : s1 + ms2 = 0, s2 â¥ 0}. Similar to Proposition 4.13 Case 2,when âr(x) = l, âf(x) = {s : s1 + ms2 â (l1 + ml2) = 0, s2 â¥ l2} andwhen âr(x) = u, âf(x) = {s : s1 + ms2 â (u1 + mu2) = 0, s2 â¥ u2}. Letâf(x) â {s : Cr(s) â¤ 0}, the case âf(x) â {s : Cr(s) â¥ 0} is analogous.Then âxâri(E)âf(x) =âxâri(E){s+âr(x) : s1 +ms2 = 0, s2 â¥ 0}= {s : l1 +ml2 â¤ s1 +ms2 â¤ u1 +mu2, Cr(s) â¤ 0}is a parabolic region.Example 4.17. For the rational function r in Example 4.8 and polytope Pwith vertices v1 = (0, 0), v2 = (1, 0) and v3 = (1, 1), let f(x) = r(x) + IP (x).From Example 4.8, we also have Cr(s) = s21 + 2s1s2 + s22 â 4s1.For an edge E = {x : x2 = x1, 0 â¤ x1 â¤ 1} between vertices v1 and v3,âr(v1) = (0, 0) and âr(v3) = (1, 1). Now âf(v1) = {s : s1 + s2 = 0, s2 â¥ 0}and âf(v2) = {s : s1 + s2â 2 = 0, s2 â¥ 1}. Also, for any x â ri(E), âf(x) â{s : Cr(s) â¥ 0}. Soâxâri(E)âf(x) = {s : 0 â¤ s1 + s2 â¤ 2, s21 + 2s1s2 + s22 â 4s1 â¥ 0}is an unbounded parabolic region. This is well illustrated in Figure 4.4 withE = E13.4.4 Structure of the conjugate domainIt is known that the conjugate of a PLQ function is a PLQ function so theimage of the domain of a PLQ function has a polyhedral subdivision [GL13].The alternate proof we provide here is useful for 2 reasons: (1) it generalizesto rational functions, and (2) it gives algorithmic insight.Theorem 4.18. For a bivariate convex quadratic function q of form (3.19)and polytope P , let f(x) = q(x) + IP (x), then the union of subdifferentialssetâxâP âf(x) has a polyhedral subdivision.594.4. Structure of the conjugate domainProof. The polytope P â R2 can be divided into vertices, faces (or edges),and the interior of P . NowâxâPâf(x) = (âªvâV âf(v))â(âªEâFâf(ri(E)))â (âªxâint(P )âf(x))(4.12)where V is the set of all vertices and F is the set of all the faces (or edges)of P . Now we compute the subdifferentials for each entity according to thefollowing three cases:Case 1(Vertices) From Lemma 4.10, for a vertex v we have âf(v) =âq(v) +NP (v), which is an unbounded polyhedral set with nonempty inte-rior.Case 2 (Edges) From Proposition 4.13, for any edge E with slope m,âxâri(E) âf(x) = {s : l1 + ml2 â¤ s1 + ms2 â¤ u1 + mu2, Cq(s) â¤ 0} is anunbounded polyhedral set with nonempty interior. In the special case whenl = u, it is a ray.Case 3 (int(P )) From Corollary 4.6,âxâint(P ) âf(x) is a contained ina line segment.Now for any v â V and E â F, âf(v)â(âxâri(E) âf(x)) is either Ïor a ray. Similarly for any E â F, (âxâri(E) âf(x))â(âxâint(P ) âf(x)) iscontained in the corresponding line segment, while for any v â V, âf(v)â(âxâint(P ) âf(x)) is a singleton {âq(v)}. SinceâxâP âf(x) = R2 is a convexset, so the subdifferential setâxâP âf(x) is the union of a finite number ofpolyhedral sets, consequently has a polyhedral subdivision.Example 4.19. For the quadratic function and the polytope in Exam-ple 4.14, we have f(x) = q(x) + IP (x). Now to computeâxâP âf(x), wecompute the subdifferential set for each entity of P .For all vertices, from Lemma 4.10, the subdifferential areâf(v1) = {s : 2s1 + 2s2 â 3 â¥ 0,â4s1 + 4s2 + 1 â¥ 0},âf(v2) = {s : 2s1 + 2s2 â 3 â¥ 0,â4s1 + 4s2 + 1 â¤ 0},âf(v3) = {s : 2s1 + 2s2 + 1 â¤ 0,â4s1 + 4s2 + 1 â¤ 0},âf(v4) = {s : 2s1 + 2s2 + 1 â¤ 0,â4s1 + 4s2 + 1 â¥ 0}.Let Eij denote the edge between vertex vi and vj . Since âf(v1) =âf(v2) and âf(v3) = âf(v4), âf(E12) and âf(E34) are contained in a ray604.4. Structure of the conjugate domainrespectively. For the remaining edgesâxâri(E23)âf(x) = {s : 1 â¤ 2s1 + 2s2 â¤ 3,â4s1 + 4s2 + 1 â¤ 0},âxâri(E14)âf(x) = {s : 1 â¤ 2s1 + 2s2 â¤ 3,â4s1 + 4s2 + 1 â¥ 0}.So the subdifferential setâxâP âf(x) is the union of all the subdifferen-tials corresponding to each edge and vertex, which are polyhedral sets. Thisis well illustrated in Figure 4.3.A simple method to obtain the polyhedral subdivision lines ofâxâP âf(x)can be formulated as follows:Step 1 For each vertex v, compute its gradient âf(v).Step 2 Connect the farthest two gradients (in terms of Euclidean distance)with a line segment. SinceâxâP âq(x) is contained in a line segment,the gradients of the rest of the vertices would lie on the segment aswell.Step 3 Compute the boundary of subdifferentials for each vertex and obtainthe subdivision lines.Corollary 4.20. For a bivariate linear function l, polytope P , let f(x) =l(x) + IP (x), then the union of subdifferentialsâxâP âf(x) has a polyhedralsubdivision.Proof. Sinceâxâint(P ) âf(x) is a singleton (from Corollary 4.5),âxâP âf(x)can be written asâvâV âf(v) which is the union of a finite number of poly-hedral sets such that for any u,w â V, âf(u) â© âf(w) is either Ï or a ray.Consequently,âxâP âf(x) has a polyhedral subdivision.Theorem 4.21. Assuming Conjecture 4.11 holds, for a bivariate rationalfunction r of form (3.18), polytope P , let f(x) = r(x) + IP (x), then theunion of the subdifferentials setâxâP âf(x) has a parabolic subdivision (seeDefinition 2.22).Proof. Let r(x) =ξ21(x)ξ2(x)+ ξ0(x) be a bivariate rational function whereξ1(x), ξ2(x) and ξ0(x) are linear functions in x. Similar to Theorem 4.18,polytope P â R2 can be divided into vertices, edges, and the interior. The614.4. Structure of the conjugate domainâf(v1)âf(v2)âf(v3)âf(v4)âxâri(E23)âf (x)âxâri(E14)âf (x)(78 , 58)(38 , 18)Figure 4.3: Polyhedral subdivision for Example 4.14624.4. Structure of the conjugate domainsetâxâP âf(x) would be the union of the subdifferentials set for each suchentity of P . It can be divided into the following three casesCase 1 (Vertices) From Lemma 4.10, for any v â V, ξ2(v) 6= 0, the setâf(v) is an unbounded polyhedral set, otherwise for some z â V such thatξ2(z) = 0, the set âf(z) is a parabolic region (from Conjecture 4.11).Case 2 (Edges) From Proposition 4.16, for any edge E â F, âxâri(E)âf(x) = {s : l1 + ml2 â¤ s1 + ms2 â¤ u1 + mu2, Cr(s) â¤ 0} is a parabolicregion. In the special case when l = u,âxâri(E) âf(x) is a ray.Case 3 (int(P )) From Corollary 4.9,âxâint(P ) âf(x) is a parabolic arc.For any v â V and E â F, âf(v)â(âxâE âf(x)) is either Ï or a ray,and (âxâri(E) âf(x))â(âxâint(P ) âf(x)) is contained in the correspondingparabolic arc. While for any v â V, âf(v)â(âxâint(P ) âf(x)) is either asingleton {âr(v)} or contained in the parabolic arc.Hence the subdifferential setâxâP âf(x) is the union of a finite numberof parabolic regions such that any two regions meet either on a parabolicarc, a ray or at a point, thus leading to a parabolic subdivision.Example 4.22. For the rational function r and polytope P in Example 4.17,f(x) = r(x) + IP (x). Now to computeâxâP âf(x), we compute the subdif-ferential set for each entity of P as follows:For vertices, from Lemma 4.10, the subdifferential at v1 and v3 areâf(v1) = {s : s1 â¤ 0, s1 + s2 â¤ 0},âf(v3) = {s : s1 + s2 â¥ 0, s2 â¥ 1}.Since v2 is the special vertex case, so âf(v2) would be computed accordingto Conjecture 4.11. However it can also be obtained by exclusion, that is,âf(v2) = R2 \u00E2xâPv2 âf(x).Similar to Example 4.22, let Eij be the edge between vertices vi andvj . Sinceâxâri(E12) âf(x) = {(0, 0)} = {âf(v1)} andâxâri(E23) âf(x) ={(1, 1)} = {âf(v3)}, âf(E12) and âf(E23) are contained in a ray respec-tively. For E13âxâri(E13)âf(x) = {s : 0 â¤ s1 + s2 â¤ 2, s21 + 2s1s2 + s22 â 4s1 â¥ 0}.So by exclusionâf(v2) ={s : s2 â¤ s1, s1 â¥ 0, s2 â¤ 1}â{s : s2 â¥ s1, s21 + 2s1s2 + s22 â 4s1 â¤ 0}.634.4. Structure of the conjugate domainâf(v1) âf(v2)âf(v3)âxâri(E13)âf(x)(1,1)(0,0)Figure 4.4: Parabolic subdivision for r and P in Example 4.22644.5. Conjugate ExpressionsThe parabolic subdivision forâxâP âf(x) is illustrated in Figure 4.4.Example 4.23. For r =36x21 + 21x1x2 + 36x22 â 81x1 + 24x2 â 252â12x1 + 9x2 + 75 andpolytope P formed by vertices v1 = (â1, 1), v2 = (â3,â3) and v3 = (â4,â3),let f(x) = r(x)+IP (x). Also from Proposition 4.7, we haveâxâdom(r) âr(x)= {s : C(s) = 0} where C(s) = 9s21 + 24s1s2 â 234s1 + 16s22 + 200s2 â 527.From Lemma 4.10, the subdifferentials at the vertices areâf(v1) = {s : 3s1 + 4s2 + 1 â¥ 0, 4s1 + 8s2 â 1 â¥ 0},âf(v2) = {s : s1 + 2s2 + 11 â¤ 0, s1 â¥ â3},âf(v3) = {s : s1 â¤ â3, 3s1 + 4s2 + 25 â¤ 0}.Again similar to Example 4.22, let Eij be an edge between vertices viand vj . Since âf(v2) = âf(v3),âxâri(E23) âf(x) would be a ray, whileâxâri(E12)âf(x) ={s : â13s1 + 4s2 â 23 â¥ 0, C(s) â¤ 0}â{s : â13s1 + 4s2â 23 â¤ 0, 4s1 + 8s2 â 1 â¤ 0, s1 + 2s2 + 11 â¥ 0},âxâri(E13)âf(x) ={s : 3s1 + 4s2 + 25 â¥ 0, 3s1 + 4s2 + 1 â¤ 0,â13s1 + 4s2â 23 â¥ 0, C(s) â¥ 0}.The parabolic subdivision for this example is shown in Figure 4.5.4.5 Conjugate Expressions4.5.1 Fractional formsLet us introduce the following fractional formsgf (s1, s2) =Ï1(s1, s2)Î¶00âÏ1/2(s1, s2)+ Ï0(s1, s2) (4.13)gq(s1, s2) = ζ11s21 + ζ12s1s2 + ζ22s22 + ζ10s1 + ζ01s2 + ζ00 (4.14)gl(s1, s2) = ζ10s1 + ζ01s2 + ζ00 (4.15)where Ï0, Ï1/2 and Ï1 are linear functions in s, and ζij â R.654.5. Conjugate Expressionsâf(v1)âf(v2)âf(v3)âxâri(E12)âf(x)âxâri(E13)âf(x)Figure 4.5: Parabolic subdivision for r and P from Example 4.5664.5. Conjugate ExpressionsTheorem 4.24. For a convex quadratic function q of form (3.19), polytopeP , let f(x) = q(x) + IP (x), then its conjugate fâ(s) has a polyhedral subdi-vision such that over each member of subdivision it has either form (4.14)or (4.15).Proof. From form (3.19), let q(x) =ξ21(x)ξ20+ ξ0(x) be a convex quadraticfunction with ξ1(x) = ξ11x1 + ξ12x2 + ξ10, ξ0(x) = ξ01x1 + ξ02x2 + ξ00 andξ20 > 0. Now the conjugate of f is computed asfâ(s) = supxâRn{ãs, xã â f(x)}.Since f(x) = q(x) + IP (x) is convex over P , so the maxima in fâ(s) canbe obtained by setting âf(x) = 0, and would lead to the following threepossible cases.Case 1 (Vertices) Let V â P be the set of all the vertices, then forv â Vfâ(s) = supxâ{v}{ãs, xã â f(x)}= ãs, xã â (q(x) + IP (x))= ãs, vã â q(v)= s1v1 + s2v2 â q(v)is a linear function in s. From Proposition 4.18, the set âf(v) is an un-bounded polyhedral set, so for any vertex v, the conjugate is a linear functiondefined on an unbounded polyhedral set.Case 2 (Edges) Let F be the set of all the edges, and E = {x : x2 =mx1 + c, l1 â¤ x1 â¤ u1} â F be an edge between vertices l and u. The casewhen m = 0 or E = {x : x1 = d, l1 â¤ x1 â¤ u1} is analogous. The conjugateexpression is then computed asfâ(s) = supxâri(E){ãs, xã â f(x)}= supxâri(E){ãs, xã â (q(x) + IP (x))}.By first order optimality condition,we havesâ (âq(x) +NP (x)) = 0.674.5. Conjugate ExpressionsFor m as the slope of the edge, NP (x) = {s : s = Î»(âm, 1), λ â¥ 0}. So wehaves1 =2ξ11ξ20t+ ξ01 âmÎ»,s2 =2ξ12ξ20t+ ξ02 + λ(4.16)where t = ξ11x1 + ξ12x2 + ξ10. Since x â ri(E), we also havex2 = mx1 + c. (4.17)From (4.16) and (4.17), we havex1 = γ1s1 + γ2s2 + γ0x2 = m(γ1s1 + γ2s2 + γ0) + cwhere all γi are defined in coefficients of q, and parameters m and c. Bysubstituting the above values in fâ(s), we getfâ(s) = ζ11s21 + ζ12s1s2 + ζ22s22 + ζ10s1 + ζ01s2 + ζ00which is a quadratic function of form (4.14) with all ζij defined in coefficientsof q(x), and parameters m and c.From Proposition 4.13, theâxâri(E) âf(x) is either an unbounded poly-hedral set with nonempty interior or a ray. So for any E, the conjugate is aquadratic function of form (4.14) defined over an unbounded polyhedral set.Whenâxâri(E) âf(x) is a ray, we do not compute the conjugate since itsvalues on the ray are given by the intersection of the conjugate expressionsof the corresponding two adjacent pieces.Case 3 (int(P )) From Corollary 4.6,âxâint(P ) âf(x) is contained in aline segment. In this case, the computation of conjugate expressions is againobtained by its neighbours.So for any convex quadratic function of form (3.19), defined over a poly-tope P , the conjugate has a polyhedral subdivision such that over eachmember of its subdivision, it is of form (4.14) or (4.15).Example 4.25. For the quadratic function q and polytope P in the Exam-ple 4.19, let f(x) = q(x) + IP (x). Similar to the domain, computation ofconjugate expressions can be divided into computing it for each entity of P .684.5. Conjugate ExpressionsFrom Theorem 4.24 case 1, the conjugate expressions for vertices arefâv1(s) =12(s1 + 2s2 â 1),fâv2(s) =316(4s1 + 4s2 â 3),fâv3(s) =116(4s1 + 4s2 â 1),fâv4(s) =s22,respectively. Now since the subdifferential set at each vertex is a polyhedralset, so for vertices, the conjugate is a linear function defined on a polyhedralset.Since from Example 4.19,âxâri(E12) âf(x) andâxâri(E34) âf(x) are rays,the conjugate expressions for the remaining edges arefâE23(s) =14(s1 + s2)2,fâE14(s) =116(4s21 + 8s1s2 + 4s22 â 4s1 + 4s2 + 1),respectively. From Example 4.19, the setsâxâri(E23) âf(x) andâxâri(E14)âf(x) are again unbounded polyhedral sets, so for the edges E23 and E14,the conjugate is a quadratic function defined on an unbounded polyhedralset.SofâP (s) =ï£±ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£³12(s1 + 2s2 â 1) s â âf(v1)316(4s1 + 4s2 â 3) s â âf(v2)14(s1 + s2)2 s â âxâri(E23) âf(x)116(4s1 + 4s2 â 1) s â âf(v3)12s2 s â âf(v4)116(4s21 + 8s1s2 + 4s22 â 4s1 + 4s2 + 1) s ââxâri(E14) âf(x),is a piecewise defined function with a polyhedral subdivision.Corollary 4.26. For a bivariate linear function l of form (3.20), a polytopeP , let f(x) = l(x) + IP (x), then its conjugate fâ(s) has a polyhedral subdi-vision such that over each member of its subdivision, it is of form (4.15).694.5. Conjugate ExpressionsProof. From Corollary 4.20,âxâP âf(x) =âvâV âf(v), so at any vertex vfâ(s) = supxâ{v}{ãs, xã â f(x)}= ãs, vã â (l(v) + IP (v))= s1v1 + s2v2 â l(v)which is a linear function in s and is of form (4.15).Theorem 4.27. Assume Conjecture 4.11 holds. For a bivariate rationalfunction r of form (3.18), a polytope P , let f(x) = r(x) + IP (x), then theconjugate fâ(s) has a parabolic subdivision such that over each member ofits subdivision it has one of the forms in (4.13)-(4.15)Proof. Let r(x) =ξ21(x)ξ2(x)+ξ0(x) be a bivariate rational function with ξ1(x) =ξ11x1+ξ12x2+ξ10, ξ2(x) = ξ21x1+ξ22x2+ξ20 and ξ0(x) = ξ01x1+ξ02x2+ξ00.The conjugate of a function f is given byfâ(s) = supxâRn{ãs, xã â f(x)}.Since f is convex (f = convP (Q)), its critical points are the maxima. Hence,similar to Theorem 4.24, computing the maxima of fâ(s) leads to the fol-lowing three cases.Case 1 (Vertices) Similar to Theorem 4.24 Case 1, by setting q = r,for any vertex vfâ(s) = s1v1 + s2v2 â r(v)is a linear function of form (4.15). In the special case, when ξ2(v) = 0,f(x) = Q(x), so fâ(s) = s1v1+s2v2âQ(v), which is again of the same form.Hence the conjugate would be a linear function of form (4.15) defined overan unbounded polyhedral set(from Lemma 4.10). In the special case, whenâf(v) is a parabolic region (Conjecture 4.11), the conjugate would again bea linear function but defined over a parabolic region.Case 2 (Edges) Similar to Theorem 4.24, let F be the set of all theedges, and E = {x : x2 = mx1 + c, l1 â¤ x1 â¤ u1} â F be an edge betweenvertices l and u, thenfâ(s) = supxâri(E){ãs, xã â f(x)}= supxâri(E){ãs, xã â (r(x) + IP (x))}.704.5. Conjugate ExpressionsBy computing the critical points, we havesâ (âr(x) +NP (x)) = 0.Again, for m as the slope of the edge, NP (x) = {s : s = Î»(âm, 1), λ â¥ 0}.So we haves1 = âÎ¾21t2 + 2ξ11t+ ξ01 âmÎ»s2 = âÎ¾22t2 + 2ξ12t+ ξ02 + λ(4.18)where t =ξ11x1 + ξ12x2 + ξ10ξ21x1 + ξ22x2 + ξ20. Since x â ri(E), we also havex2 = mx1 + c. (4.19)From (4.18) and (4.19), we havex1 =ï£±ï£´ï£²ï£´ï£³Î³10s1 + γ01s2 + γ00 when ξ21 +mξ22 = 0γ00 Â± Î³1/2âÎ³10/2s1 + γ01/2s2 + Î³00/2±Î³â1/2âÎ³10/2s1 + γ01/2s2 + γ00/2otherwise,(4.20)where all γij and γij/k are defined in the coefficients of r, and parameters mand c. When ξ21 +mξ22 6= 0, solving (4.18) and (4.19), leads to a quadraticequation in t with coefficients as linear functions in s, and so Â± sign comesas a result of the quadratic roots.By subtituting (4.20) and (4.19) in fâ(s), when ξ21 +mξ22 6= 0, we havefâ(s) =Ï1(s1, s2)Î¶00âÏ1/2(s1, s2)+ Ï0(s1, s2),and when ξ21 +mξ22 = 0,fâ(s) = ζ11s21 + ζ12s1s2 + ζ22s22 + ζ10s1 + ζ01s2 + ζ00,where all ζij , Ïi and Ïi/j are defined in the coefficients of r, and paramatersm and c. Also, Ïi(s) and Ïi/j(s) are linear functions in s. A detailedformulation is provided in Appendix C.From Proposition 4.16,âxâri(E) âf(x) is either a parabolic region or aray. So for any E, the conjugate is a fractional function of form (4.13) definedover a parabolic region. Similar to Theorem 4.24 case 2, whenâxâri(E) âf(x)is a ray, the computation of the conjugate is deduced from its neighboursby continuity.714.5. Conjugate ExpressionsCase 3 (int(P )) Similar to Theorem 4.24, sinceâxâint(P ) âf(x) is con-tained in a parabolic arc (from Corollary 4.9), so the computation of theconjugate is deduced by continuity.So for any rational function of form (3.19), defined over a polytope P ,its conjugate has a parabolic subdivision such that over each member of thesubdivision it has one of the forms in (4.13)- (4.15).Example 4.28. For the rational function r and polytope P in Example 4.23,f(x) = r(x) + IP (x). Similar to Example 4.25, we obtain the full conjugateby computing the conjugate for the vertices and edges of P only.For the vertices, from Theorem 4.27 case 1, the conjugate expressionsfor vertices arefâv1(s) = âs1 + s2 + 1,fâv2(s) = â3s1 â 3s2 â 9,fâv3(s) = â4s1 â 3s2 â 12,and from Example 4.23, the subdifferential set for each vertex is an un-bounded polyhedral set. So for each vertex of P , the conjugate is a linearfunction defined over an unbounded polyhedral set.Asâxâri(E23) âf(x) is a ray, for the remaining edgesfâE12(s) =336s1 + 672s2 â 12432â3(âs1 â 2s2 + 37)â 17s1 â 31s2 + 1181,fâE13(s) =148(9s21 + 24s1s2 + 16s22 â 42s1 + 56s2 + 49)respectively. So from Example 4.23, for edge E12, the conjugate is a frac-tional function of form (4.13) with parabolic region, while for edge E13, theconjugate is a quadratic function defined over a parabolic region again. SofâP (s) =ï£±ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£³âs1 + s2 + 1 s â âf(v1)336s1+672s2â12432â3(âs1â2s2+37)â 17s1 â 31s2 + 1181 s ââxâri(E12) âf(x)â3s1 â 3s2 â 9 s â âf(v2)â4s1 â 3s2 â 12 s â âf(v3)(9s21+24s1s2+16s22â42s1+56s2+49)48 s ââxâri(E13) âf(x)is a piecewise defined function with parabolic subdivision.72Chapter 5Algorithmic computation ofthe Conjugate for a class ofBivariate Rational functionsover a polytopeThis section covers the data structures and the algorithmic design usedfor the computation of the conjugate. We start with our proposed divisionof a parabolic region and a hybrid data structure between numeric andsymbolic forms to store the domain and the expressions of the conjugate.After that, we discuss our algorithmic design along with two examples forthe computation of the conjugate.5.1 Data StructuresThe input data structure is similar to the output data structure in Sec-tion 3.6.1. The coefficients of the rational function are stored in an array andthe regions in an N Ã 2 vertex or N Ã 4 halfspaces matrix. The output datastructure chosen for this module is a hybrid between numerical and symbolicstructures. Since the conjugate of a rational function is a piecewise fractionalfunction with a parabolically divided domain, the fractional expressions arestored in a symbolic array and the domain as a cell of matrices.Definition 5.1. Cell (Data Structure). A cell is tabular form of data struc-ture with rows and columns based indexing for its elements where each ele-ment can again be a cell, matrix or vector. The advantage of a cell is thatelements can be of different sizes and can have a different data type.Example 5.2. A cell C of size 3Ã 2 could beC{1, 1} = [1, 3, 2, 4] C{1, 2} = [x1 + x2 + 1]735.1. Data StructuresFigure 5.1: A subdifferential regionC{2, 1} =[2x1 + 3x2 + 5 = 0s21 + 3s1s2 + 7s22 + 9s2 = 0]C{2, 2} = 227C{3, 1} =[0.1 0.70.5 .47]C{3, 2} = [ ]where x1, x2, s1 and s2 are symbolic variables.5.1.1 Parabolic regionTo store the conjugate domain for a rational function, we need a wayto store the subdifferentials with parabolic regions. Unlike the polyhedralcase, storing all the active inequalities does not lead to the same region.For example, in Figure 5.1, storing only the active inequalities leads to theparabolic region with only one active inequality as shown in Figure 5.2.To overcome this problem, we store the region with a parabolic arc asa parabolic region with parabolic inequality and a halfspace obtained byconnecting the endpoints of the arc, and the rest of the region as a polyhedralregion. This is well illustrated in Figure 5.3.745.1. Data StructuresFigure 5.2: Region with only active constraintsFigure 5.3: Proposed subdivision for storing a subdifferential region755.1. Data StructuresTable 5.1: Conjugate domain for a rational functionEntity Form 1 Form 2Vertex Polyhedral Parabolic + halfspaceâPolyhedralri(E) Parabolic + Polyhedral Parabolic + halfspaceâPolyhedralint(P ) Parabolic arc â5.1.2 Output Data StructureFor the conjugate with N pieces, the domain is stored as an N Ã 2cell. The first column contains all the halfspaces forming a region while thesecond column contains all the parabolic regions in the symbolic form. Foreach row, the region is given by the intersection of polyhedral and parabolicregions in first and second columns. Each region can have up to a maximumof only one parabolic region. The possible structures for the domain of theconjugate of a rational function of form (3.18) is shown in Table 5.1.2. Forfâ(s) =ï£±ï£´ï£´ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£´ï£´ï£³f1(s) {s : A1s â¤ b1, Cr(s) â¤ 0}f2(s) {s : A2s â¤ b2}... ...... ...fN (s) {s : ANs â¤ bN , Cr(s) â¥ 0},the conjugate domain is stored ascdiv{1, 1} =ï£®ï£°1 m11 c11 11 m12 c12 0. . . .ï£¹ï£» cdiv{1, 2} = [Cr(s) â¤ 0]cdiv{2, 1} =ï£®ï£°1 m21 c21 11 m22 c22 0. . . .ï£¹ï£» cdiv{2, 2} = []· · · · · Â·cdiv{N, 1} =ï£®ï£°1 mN1 cN1 11 mN2 cN2 0. . . .ï£¹ï£» cdiv{N, 2} = [Cr(s) â¥ 0],while the expressions are stored in a symbolic array ascxps =[f1(s) f2(s) · · · fN (s)].This is well illustrated in Figure 5.4.765.1. Data Structurescdiv[ ][ ][ ][ Cr â¤ 0 ][ Cr â¥ 0 ]P1P2P3P4P5P1P2P3P4P5Crcxps f1(s) f2(s) â¯ â¯ fN(s)Figure 5.4: Output Data structure775.1. Data StructuresExample 5.3. For a rational function r(x) =x2x2 â x1 + 1 and a polytope Pwith vertices v1 = (0, 0), v2 = (1, 0) and v3 = (1, 1), let f(x) = r(x) + IP (x).Then the conjugate of f isfâP (s) =ï£±ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£³0 s â âf(v1)s1 s â âf(v2)s1 + s2 â 1 s â âf(v3)14(s1 + s2)2 s â âxâri(E13) âf(x)where the subdifferential sets are defined in Example 4.22. So it is stored ascxps =[0 s1 s1 s1 + s2 â 1 14(s1 + s2)2],andcdiv{1, 1} =[1 â1 2 11 0 1 1]cdiv{1, 2} = [ ]cdiv{2, 1} = [1 1 0 1] cdiv{2, 2} = [ s21 + 2s1s2 + s22 â 4s1 â¤ 0 ]cdiv{3, 1} =ï£®ï£°1 1 0 01 0 1 00 1 0 0ï£¹ï£» cdiv{3, 2} = [ ]cdiv{4, 1} =[0 1 0 11 â1 0 0]cdiv{4, 2} = [ ]cdiv{5, 1} =ï£®ï£°1 â1 0 11 â1 2 01 1 0 1ï£¹ï£» cdiv{5, 2} = [ s21 + 2s1s2 + s22 â 4s1 â¥ 0 ].Since the domain of the conjugate of a quadratic or linear function onlyhas a polyhedral subdivision as shown in Table 5.1.2, their domains with Npieces are stored as an N Ã 1 cell. The symbolic column is not required inthis case.Example 5.4. For a quadratic function q(x) =14(x1 + x2)2 +18(x1 â x2)and polytope P with vertices v1 = (0.5, 1), v2 = (0.75, 0.75), v3 = (0.25, 0.25)785.2. Algorithmand v4 = (0, 0.5), let f(x) = q(x) + IP (x). ThenfâP (s) =ï£±ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£³12(s1 + 2s2 â 1) s â âf(v1)316(4s1 + 4s2 â 3) s â âf(v2)14(s1 + s2)2 s â âxâri(E23) âf(x)116(4s1 + 4s2 â 1) s â âf(v3)12s2 s â âf(v4)116(4s21 + 8s1s2 + 4s22 â 4s1 + 4s2 + 1) s ââxâri(E14) âf(x),where the subdifferential sets are defined in Example 4.19. It is stored ascxps = [s1/2 + s2 â 1/2, (3s1)/4 + (3s2)/4â 9/16, (s1 + s2)2/4,s1/4 + s2/4â 1/16, s2/2, s21/4 + (s1s2)/2â s1/4 + s22/4 + s2/4 + 1/16]andcdiv{1, 1} =[1 â1 1.50 11 1 â0.25 1]cdiv{2, 1} =[1 1 â0.25 01 â1 1.50 1]cdiv{3, 1} =ï£®ï£°1 â1 1.50 01 â1 0.50 11 1 â0.25 0ï£¹ï£»cdiv{4, 1} =[1 â1 0.50 01 1 â0.25 0]cdiv{5, 1} =[1 1 â0.25 11 â1 0.50 0]cdiv{6, 1} =ï£®ï£°1 â1 0.50 11 â1 1.50 01 1 â0.25 1ï£¹ï£» .5.2 AlgorithmThis section covers the main algorithms in the module to compute theconjugate of bivariate functions of forms (3.18)-(3.20). The main algorithmfollows the same method as Theorem 4.27. We also discuss the algorithmused to compute the subdifferentials for the edges of the polytope P .795.2. AlgorithmTable 5.2: Conjugate domain for quadratic functionsEntity FormVertex Polyhedralri(E) Polyhedralint(P ) Line SegmentAlgorithm 6 Computing the conjugate for rational function1: function compute conjugate(rx,P)2: [cdiv, cdivt,ncrgn] = compute conjugate domain div(rx, P);3: cxpst = compute conjugatexps vertices(ncrgn);4: cxps = compute conjugatexps edges(cdivt,cxpst);5: return cxps,cdiv6: end function5.2.1 Main AlgorithmThe main algorithm is presented in Algorithm 6. It takes the rationalfunction rx and the polytope region P following the input data structures asinput arguments. We start with computing the domain division and store itin the variable cdiv. The variables cdivt and ncrgn contain the intermediatedata for computing the domain division and the normal cone regions respec-tively. The variable cdivt follows the same output data structure, however,it stores the mapping between primal polytope and the parabolic or polyhe-dral subdivision in the dual domain. The variable ncrgn stores the verticesand their normal cones in the primal as an N Ã 2 cell where N is the sumof the total number of vertices and edges of polytope P.By using ncrgn, we compute the conjugate expressions corresponding tothe vertices of P and store it in the variable cxpst. After that, we computethe remaining expressions using the variables cdivt and cxpst, and store itin cxps. The variables cxps and cdiv contains the conjugate expressions anddomain division with the relation being mapped as the indices of the arraycxps to the rows in cdiv.The subroutine compute conjugate domain div() has been presented inAlgorithm 7. It iterates through each entity once and computes their sub-differentials using the different subroutines. Since it traverse through themonly once, for a polytope with n vertices, the time complexity is T (n) =c1 + c2n + c3n + c4n where ci are constants, and so is bounded by O (n).Similarly, for computing the expressions, we go through each entity again805.2. AlgorithmAlgorithm 7 Computing the domain of the conjugate1: function compute conjugate domain div(rx, rgn)2: rels1s2 = compute curve equation(rx,rgn); . equation of the curve3: ncrgn = compute normal cone for polyhedral domain(rgn);4: cdivt = compute subdiff vertices(ncrgn);5: [cdiv,iv] = compute subdiff edges(rels1s2, cdivt, ncrgn);6: return cdiv, cdivt, ncrgn, iv7: end functiononly once, then overall the time complexity of the algorithm remains O (n).5.2.2 Algorithm for edgesAn algorithm for computing the union for subdifferentials at the edgeshas been shown in Algorithm 8. For each edge, we start with finding itscorresponding adjacent vertices and the parallel halfspaces forming the sub-differentials region. We compute the normal vector on the edge and storeit in variable nv. After that, using the equation of the curve rels1s2, theadjacent vertices, parallel halfspaces prs and the normal vector, we computethe parabolic region with the ââ¤ 0â inequality. To verify that this is ourparabolic region, we compute the gradient at the midpoint of the edge andthen add the normal vector to it. Now if this point satisfies all the inequali-ties of the ââ¤ 0â parabolic region, we set the flag islessthan to true else it isfalse. When flag is true, it leads to a subdivision of the region to a parabolicand polyhedral regions, otherwise the parabolic region is ââ¥ 0â and can berepresented with only one parabolic region.Similarly while computing the expressions for the union of subdifferen-tials on the edges, the relation between x and s has a quadratic form andleads to two relations as shown in 4.27 [Case 2 edges]. Now since the conju-gate would be continuous on the boundary with neighboring expressions, wefind the conjugate expression by substituting the equation of the boundariesin the two and verifying its equality with the neighbors. The case of thequadratic and linear functions are simple and do not require any furthersubdivision or verification of the expressions.815.2. AlgorithmAlgorithm 8 Computing the subdifferentials for the edges1: function compute subdiff edges(rels1s2, cdivt, ncrgn)2: for i = 1 to size(cdivt,1) do3: if isempty(cdivt{i,1}) then . if an edge4: [vli,vri] = choose proper indexpm1(i,size(ncrgn,1));5: v = [cdivt{vli,1}; cdivt{vri,1}];6: hslv = cdivt{vli,2}(2,:); hsrv = cdivt{vri,2}(1,:);7: prs = [hslv(1:3),abs(1-hslv(4)); hsrv(1:3),abs(1-hsrv(4))];8: [nv] = compute normal vector(ncrgn{i,1},ncrgn{vli,1},ncrgn{i,2}(1,:));9: [ldr] = compute curve less than region parll(rels1s2, v, prs,nv);10: [x1v,x2v] = compute midpoint on edge(ncrgn{i,1},ncrgn{vli,1}, ncrgn{vri,1});11: sc = [drx1(x1v,x2v), drx2(x1v,x2v)];12: ptc = sc + nv;13: if inside polytope check(ldr{1,1},ptc(1),ptc(2)) andrl(ptc(1),ptc(2))â¤ 0 then14: islessthan = true;15: else if inside polytope check(ldr{2,1},ptc(1),ptc(2)) then16: islessthan = true;17: else18: islessthan = false;19: end if20: if islessthan then21: cdiv{j,1} = ldr{1,1}; cdiv{j,2} = ldr{1,2}; iv(j) = i; j++;22: cdiv{j,1} = ldr{2,1}; iv(j) = i; j++;23: else24: cdiv{j,1} = [prs; ldr{1,1}]; cdiv{j,2} = [rels1s2 >= 0];25: iv(j) = i; j++;26: end if27: end if28: end for29: return cdiv, iv30: end function825.3. Example 1Figure 5.5: Normal cone division lines for Example 5.35.3 Example 1For a bivariate rational function r(x) =x22x2 â x1 + 1 defined over a poly-tope P with vertices v1 = (1, 1), v2 = (1, 0) and v3 = (0, 0), let f(x) =r(x) + IP (x). We start with computing the normal cone division lines forthe polytope as shown in Figure 5.5 and store it in variable ncrgn asncrgn{1, 1} = [1 1] ncrgn{1, 2} = [1 â1 2 11 0 1 1]ncrgn{2, 1} = [0 1 â1 1] ncrgn{2, 2} = [1 0 1 01 0 0 1]ncrgn{3, 1} = [1 0] ncrgn{3, 2} = [1 0 0 00 1 â1 0]ncrgn{4, 1} = [1 0 0 1] ncrgn{4, 2} = [0 1 â1 10 1 0 0]ncrgn{5, 1} = [0 0] ncrgn{5, 2} = [0 1 0 11 â1 0 0]ncrgn{6, 1} = [1 1 0 0] ncrgn{6, 2} = [1 â1 0 11 â1 2 0].All the arrays in the first column represent an entity corresponding to835.3. Example 1the polytope P with adjacent elements being the neighboring entity whentraversed clockwise direction while the second row shows the normal cone(or union of normal cones) region corresponding to that entity. The arrays oflength 2 are vertices while ones with length 4 are the equation of the edges.For example, the equation of line given by ncrgn{2, 1} is the relation ofthe points lying on the edge between vertices ncrgn{1, 1} and ncrgn{3, 1},while the regions given by ncrgn{1, 2} is the normal cone for vertex v3 andncrgn{2, 2} is the union of normal cones when presented in primal domain.We compute the gradients of r asdrx1(x1, x2) =(x2x2 â x1 + 1)2drx2(x1, x2) = 2(x2x2 â x1 + 1)â(x2x2 â x1 + 1)2where drx1 =âr(x)âx1and drx2 =âr(x)âx2are symbolic variables storing thegradients respectively. Using the gradients, we compute the equation of theparabolic curve for r, as shown in Example 4.8, and store it in the symbolicvariable rels1s2 = s21 + 2s1s2 + s22 â 4s1.After that, we compute the value of gradients at all the vertices as theendpoints of the parabolic arc would be the farthest two vertices. Nowdrx1(v1) = 1 drx2(v1) = 1,drx1(v3) = 0 drx2(v3) = 0,however, since r(v2) =00is not defined, andlimx1âv21(limx2âv22âr(x)âx1)6= limx2âv22(limx1âv21âr(x)âx1)andlimx1âv21(limx2âv22âr(x)âx2)6= limx2âv22(limx1âv21âr(x)âx2),the set âf(v2) is a union of parabolic and polyhedral region given by R2 \u00E2xâPv2 âf(x), so would be computed at the end. Also for edges E12 andE23,drx1(E12) = drx1(v1) and drx2(E12) = drx2(v1),drx1(E23) = drx1(v3) and drx2(E23) = drx2(v3),845.3. Example 1the subdifferentials for these edges lie along a ray and can be deduced by theneighboring entities. Along with computing the gradients and the normalcones at vertices, we set up an identifier for v2 in the transitioning variablecdivt ascdivt{1, 1} = [1 1] cdivt{1, 2} = [1 1 2 11 0 1 1]cdivt{2, 1} = [NaN NaN NaN NaN] cdivt{2, 2} = [ ]cdivt{3, 1} = [NaN NaN] cdivt{3, 2} = [1 0 1 00 1 0 0]cdivt{4, 1} = [NaN NaN NaN NaN] cdivt{4, 2} = [ ]cdivt{5, 1} = [0 0] cdivt{5, 2} = [0 1 0 11 â1 0 0]cdivt{6, 1} = [ ] cdivt{6, 2} = [ ].The length of arrays in the first column with two NaN values signifies thevertex with the special case and the one with four NaN values signifies anedge for which computation of subdifferentials is not required as it is a rayand can be deduced by the neighboring entities.At this point, we only need to compute âf(v2) andâxâri(E13) âf(x).For the edge, we start with finding the gradients at the adjacent vertices,which in this case are âf(v1) = [0, 0] and âf(v3) = [1, 1]. Now the parallelhalfspaces forming the subdifferential regions would be the opposite regionof one of the halfspaces forming the normal cones at these vertices, that is,prs =[1 â1 0 11 â1 2 0]which is [ cdivt{5, 1}(2, 1 : 3), abs(1 â cdivt{5, 1}(2, 4)); cdivt{1, 1}(1, 1 :3), abs(1 â cdivt{1, 1}(1, 4)) ]. After that, we compute the normal vectorcorresponding to the edge of P and store it as nv = [â1/â2, 1/â2]. Thenormal cone at any point x â ri(E13) is Î»Ãnv for λ â¥ 0, so s = âf(x)+nv ââxâri(E13) âf(x). We pick x to be the mid point of the vertices v1 and v3lying on the edge and store it in mp = [1/2, 1/2]. Since we do not knowwhich side of the parbolic arc the region lies, we assume that the regionlies on the ââ¤ 0â side of the parabolic arc, which we verify by the criteria855.3. Example 1s â âxâri(E13) âf(x). This region is stored asldr{1, 1} = [1 1 0 1] ldr{1, 2} = [ s21 + 2s1s2 + s22 â 4s1 â¤ 0 ]ldr{2, 1} =ï£®ï£°1 â1 0 11 â1 2 01 1 0 0ï£¹ï£» ldr{2, 2} = [ ].However, since s = âf(mp) + nv â [â0.4571, 1.4571 ] fails to satisfy theconstraints for ldr ass21 + 2s1s2 + s22 â 4s1 = 2.8284  0,it belongs to the ââ¥ 0â side of the parabolic arc, so the subdifferentials forthe region is stored ascdiv{5, 1} =ï£®ï£°1 â1 0 11 â1 2 01 1 0 1ï£¹ï£» cdiv{5, 2} = [ 0 â¤ s21 + 2s1s2 â 4s1 + s22 ].Now at the end, we compute âf(v2). This is identified by an array withtwo NaN values in cdivt. Since âf(v2) is a convex set, it belongs to theââ¤ 0â side of the parabolic, we only need to compute the halfspaces formingthe boundary of this region. Since âf(E12) and âf(E23) maps to âf(v1)and âf(v3) respectively, similar to the edge E13, the halfspaces would bethe opposite side of the sharing hyperplane from the normal cones at theadjacent vertices, and would be stored as a subdivision with one parabolicand one polyhedral region ascdiv{2, 1} = [1 1 0 1] cdiv{2, 2} = [ s21 + 2s1s2 + s22 â 4s1 â¤ 0 ]cdiv{3, 1} =ï£®ï£°1 1 0 01 0 1 00 1 0 0ï£¹ï£» cdiv{3, 2} = [ ]Finally, all the regions would be stored in cdiv ascdiv{1, 1} =[1 â1 2 11 0 1 1]cdiv{1, 2} = [ ]cdiv{2, 1} = [1 1 0 1] cdiv{2, 2} = [ s21 + 2s1s2 â 4s1 + s22 â¤ 0 ]cdiv{3, 1} =ï£®ï£°1 1 0 01 0 1 00 1 â0 0ï£¹ï£» cdiv{3, 2} = [ ]865.3. Example 1cdiv{4, 1} =[0 1 â0 11 â1 0 0]cdiv{4, 2} = [ ]cdiv{5, 1} =ï£®ï£°1 â1 0 11 â1 2 01 1 0 1ï£¹ï£» cdiv{5, 2} = [ 0 â¤ s21 + 2s1s2 â 4s1 + s22 ].Along with storing the transition matrix, we also store the indices mappingto compute the corresponding expression for the subregions in cdiv as anarrayiv = [1, 3, 3, 5, 6]where the index of iv corresponds to the row in cdiv and the value at theindex to the row in cdivt.Next we use the variable ncrgn and rx to compute the conjugate expres-sions for vertices v1 and v3 ascxpst(1) = s1 + s2 â 1,andcxpst(5) = 0.For v2, f(v2) =00is undefined, however, since f is the closed convex en-velope of a nonconvex quadratic function Q, limxâv2f(x) should exist, henceupon computing limx1âv21(limx2âv22f(x))= 0. So the conjugate expressioncorresponding to v2 iscxpst(3) = 1 · s1 + 0 · s2 â 0 = s1.Now for the edge E12, by computing the critical point, from Theo-rem 4.27 [Case 2], we have sâ (âr(x) +NP (x)) = 0, sos1 = t2 + λs2 = 2tâ t2 â λwhere t =x2x2 â x1 + 1. By adding both of the above equations, we gett =s1 + s22.875.3. Example 1Figure 5.6: Conjugate for Example 5.3Since the points lie on the edge E13, we also have x2 = x1. So x2 = x1 =s1 + s22. By using the optimal values of x1 and x2, upon simplification, theconjugate expression becomescxpst(6) =(s1 + s2)24.Now by using cxpst and iv, we get cxps = cxpst(iv) ascxps =[s1 + s2 â 1, s1, s1, 0, (s1 + s2)24]The conjugate in mathematical form isfâP (s) =ï£±ï£´ï£´ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£´ï£´ï£³0 s â R1s1 s â R2s1 s â R3s1 + s2 â 1 s â R414(s1 + s2)2 s â R5885.4. Example 2whereR1 = {s : s2 â¥ âs1 + 2, s2 â¥ 1}R2 = {s : s2 â¥ s1, s21 + 2s1s2 â 4s1 + s22 â¤ 0}R3 = {s : s2 â¤ s1, s2 â¤ 1, s1 â¥ 0}R4 = {s : 0 â¤ s1, s2 â¤ âs1}R5 = {s : s2 â¥ âs1, s2 â¤ âs1 + 2, s2 â¥ s1, s21 + 2s1s2 â 4s1 + s22 â¥ 0}.It is shown in Figure 5.6.5.4 Example 2For r(x) =36x21 + 21x1x2 + 36x22 â 81x1 + 24x2 â 252â12x1 + 9x2 + 75 and polytope Pformed by vertices v1 = (â1, 1), v2 = (â3,â3) and v3 = (â4,â3), let f(x) =r(x)+IP (x). We start with computing the normal cones for each entity andstore them in variable ncrgn asncrgn{1, 1} = [â1 1] ncrgn{1, 2} = [1 â0.75 0.25 11 â0.50 0.50 1]ncrgn{2, 1} = [1 2 3 1] ncrgn{2, 2} = [1 â0.50 0.50 01 â0.50 â4.50 1]ncrgn{3, 1} = [â3 â3] ncrgn{3, 2} = [1 â0.50 â4.50 00 1 3 0]ncrgn{4, 1} = [1 0 â3 1] ncrgn{4, 2} = [0 1 3 10 1 4 0]ncrgn{5, 1} = [â4 â3] ncrgn{5, 2} = [0 1 4 11 â0.75 â6 0]ncrgn{6, 1} = [1 1.33 2.33 0] ncrgn{6, 2} = [1 â0.75 â6 11 â0.75 0.25 0]The normal cone division lines is shown in Figure 5.7. Next we compute thegradients asdrx1(x) =â48x21 + 72x1x2 + 600x1 + 69x22 + 126x2 â 1011(3x2 â 4x1 + 25)2drx2(x) =â64x21 â 96x1x2 + 224x1 + 36x22 + 600x2 + 452(3x2 â 4x1 + 25)2 .895.4. Example 2Figure 5.7: Normal cone division lines for Example 5.4Using the above gradients, we compute the equation of the parabolic curvefor r and store it in the symbolic variable rl(s1, s2) = 9s21 +24s1s2â234s1 +16s22 + 200s2 â 527.Now the gradients at each vertex isdrx1(v1) = â1.5 drx2(v1) = 0.8750drx1(v2) = â3 drx2(v2) = â4drx1(v3) = â3 drx2(v3) = â4.Since the âf(v2) = âf(v3), the set âf(E23) is a ray. So the transitionvariable cdivt iscdivt{1, 1} = [â1.50 0.88] cdivt{1, 2} = [1 â0.75 â0.25 11 â0.50 0.12 1]cdivt{2, 1} = [ ] cdivt{2, 2} = [ ]cdivt{3, 1} = [â3 â4] cdivt{3, 2} = [1 â0.50 â5.50 00 1 3 0]cdivt{4, 1} = [nan, nan, nan, nan] cdivt{4, 2} = [ ]cdivt{5, 1} = [â3 â4] cdivt{5, 2} = [0 1 3 11 â0.75 â6.25 0]cdivt{6, 1} = [ ] cdivt{6, 2} = [ ]905.4. Example 2where the nan values in cdivt{4, 1} represents the subdifferentials infor-mation E23 which can be deduced by the neighboring entities, so directcomputation of the subdifferentials is not required.At this point, we only need to computeâxâri(E12) âf(x) andâxâri(E13)âf(x). For edge E12, we compute the parallel halfspaces forming the sub-differentials region, which are the regions on the opposite side of one thehalfspaces forming the normal cones at the adjacent vertices. We store thisin variableprs =[1 â0.50 0.12 01 â0.50 â5.50 1].Now we only need to verify on which side of the parabolic arc does thesubdifferentials region lie. So we formulate the ââ¤ 0â side of the inequalityregion asldr{1, 1} = [1 3.25 5.75 1] ldr{1, 2} = [ 9s21 + 24s1s2 â 234s1+ 16s22 + 200s2 â 527 â¤ 0 ]ldr{2, 1} =ï£®ï£°1 â0.50 0.12 01 â0.50 â5.50 11 3.25 5.75 0ï£¹ï£» ldr{2, 2} = [ ].Then we compute the normal vector nv = [ 0.8944, â0.4472 ] and thevalue of the gradient at the mid point of the vertices and on the edge assc = [ â2.5733, â1.1200 ]. Next, we verify if the point s = nv + sc =[ â1.6789, â1.5672 ] lie in the region ldr. So9s21 + 24s1s2 â 234s1 + 16s22 + 200s2 â 527 = â319.7627 â¤ 0,hence the subdifferential region corresponding to the edge E12 is ldr and itwould be stored in cdiv ascdiv{2, 1} = [1 3.25 5.75 1] cdiv{2, 2} = [ 9s21 + 24s1s2 â 234s1+ 16s22 + 200s2 â 527 â¤ 0 ]cdiv{3, 1} =ï£®ï£°1 â0.50 0.12 01 â0.50 â5.50 11 3.25 5.75 0ï£¹ï£» cdiv{3, 2} = [ ].Similarly for the edge E13, the parallel halfspaces would beprs =[1 â0.75 â6.25 11 â0.75 â0.25 0],915.4. Example 2and nv = [ â0.8, 0.6 ] and sc = [ â2.6250, â1.2812 ]. So s = [ â3.4250,â0.6813 ], but9s21 + 24s1s2 â 234s1 + 16s22 + 200s2 â 527 = 307.2  0,so it lies on the ââ¥ 0â side of the inequality. The subdifferential region, inthis case, is then stored ascdiv{6, 1} =ï£®ï£°1 â0.75 â6.25 11 â0.75 â0.25 01 3.25 5.75 1ï£¹ï£»cdiv{6, 2} = [ 9s21 + 24s1s2 â 234s1 + 16s22 + 200s2 â 527 â¥ 0 ].Finally by all the subdifferential regions are stored in cdiv ascdiv{1, 1} =[1 â0.75 â0.25 11 â0.50 0.12 1]cdiv{1, 2} = [ ]cdiv{2, 1} = [1 3.25 5.75 1] cdiv{2, 2} = [ 9s21 + 24s1s2 â 234s1+ 16s22 + 200s2 â 527 â¤ 0 ]cdiv{3, 1} =ï£®ï£°1 â0.50 0.12 01 â0.50 â5.50 11 3.25 5.75 0ï£¹ï£» cdiv{3, 2} = [ ]cdiv{4, 1} =[1 â0.50 â5.50 00 1 3 0]cdiv{4, 2} = [ ]cdiv{5, 1} =[0 1 3 11 â0.75 â6.25 0]cdiv{5, 2} = [ ]cdiv{6, 1} =ï£®ï£°1 â0.75 â6.25 11 â0.75 â0.25 01 3.25 5.75 1ï£¹ï£» cdiv{6, 2} = [ 9s21 + 24s1s2 â 234s1+ 16s22 + 200s2 â 527 â¥ 0 ]and the index mapping variable asiv = [1, 2, 2, 3, 5, 6].Now using ncrgn and r, the conjugate expressions corresponding to thevertices arecxpst(1) = âs1 + s2 + 1,cxpst(3) = â3s1 â 3s2 â 9,andcxpst(5) = â4s1 â 3s2 â 12.925.4. Example 2For the edge E12, by computing the critical point, from Theorem 4.27[Case 2], we have sâ (âr(x) +NP (x)) = 0, sos1 = drx1(x) + α (5.1)s2 = drx2(x)â α/2 (5.2)and since the points lie on an edge, we also havex2 = 2x1 + 3. (5.3)By eliminating α from the Equations (5.1) - (5.3), we have(s1 + 2s2 â 37)x21 + (34s1 + 68s2 â 1258)x1 + 289s1 + 578s2 â 1285 = 0,and upon simplification, we getx1 =â17s1 â 34s2 + 629± 56â3 · (âs1 â 2s2 + 37)s1 + 2s2 â 37 .Now by using the above equation and (5.3), we have two conjugate expres-sions astcxp(1) =336s1 + 672s2 â 12432â3(âs1 â 2s2 + 37)+ 17s1 + 31s2 â 1181tcxp(2) =336s1 + 672s2 â 12432â3(âs1 â 2s2 + 37)â 17s1 â 31s2 + 1181.Since the conjugate is continuous, we verify the continuity of the expressionsat the boundaries with the neighboring entites. Now the neighbouring enti-ties in this case are the subdiffential regions corresponding to vertices v1 andv2, and the shared boundaries are the rays s2 = â0.5s1 + 0.1250, s1 â¥ â1.5and s2 = â0.5s1 â 5.5, s1 â¥ â3 with expressions cl(s) = âs1 + s2 + 1and cr(s) = â3s1 â 3s2 â 9 respectively. So on boundaries, clb(s1) =cl(s1,â0.5s1 + 0.1250) = â32s1 + 98 and crb(s1) = cr(s1,â0.5s1 â 5.5) =â32 s1 +152 . Now the restrictions of tcxp(1) along the boundaries aretpcl(1) =â32s1 +188258, tpcr(1) =â32s1 +53912andtpcl(2) = â32s1 +98, tpcr(2) =â32s1 +152.935.4. Example 2So for edge E12, the conjugate expression iscxpst(2) =336s1 + 672s2 â 12432â3(âs1 â 2s2 + 37)â 17s1 â 31s2 + 1181.Similarly for edge E13, for sâ (âr(x) +NP (x)) = 0, we haves1 = drx1(x) + α (5.4)s2 = drx2(x)â 3α4. (5.5)Also from the equation of the line forming the edgex2 =4x13+73. (5.6)Now by eliminating α from the Equations (5.4) - (5.5), we getx1 =38s1 +s22â 78.Finally, by using the above equation and (5.6), the conjugate expression forE13 comes out to becxpst(6) =148(9s21 + 24s1s2 + 16s22 â 42s1 + 56s2 + 49).So by using cxpst and iv, the conjugate expressions cxps = cxpst(iv) arecxps =[â s1 + s2 + 1, 336s1 + 672s2 â 12432â3(âs1 â 2s2 + 37)â 17s1 â 31s2 + 1181,336s1 + 672s2 â 12432â3(âs1 â 2s2 + 37)â 17s1 â 31s2 + 1181, â3s1 â 3s2 â 9,â 4s1 â 3s2 â 12, 148(9s21 + 24s1s2 + 16s22 â 42s1 + 56s2 + 49)].The conjugate in mathematical form isfâP (s) =ï£±ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£´ï£³âs1 + s2 + 1 s â R1336s1 + 672s2 â 12432â3(âs1 â 2s2 + 37)â 17s1 â 31s2 + 1181 s â R2336s1 + 672s2 â 12432â3(âs1 â 2s2 + 37)â 17s1 â 31s2 + 1181 s â R3â3s1 â 3s2 â 9 s â R4â4s1 â 3s2 â 12 s â R5148(9s21 + 24s1s2 + 16s22 â 42s1 + 56s2 + 49) s â R6945.4. Example 2Figure 5.8: Conjugate for Example 5.4whereR1 = {s : 3s1 + 4s2 + 1 â¥ 0, 4s1 + 8s2 â 1 â¥ 0}R2 = {s : â13s1 + 4s2 â 23 â¥ 0, C(s) â¤ 0}R3 = {s : â13s1 + 4s2 â 23 â¤ 0, 4s1 + 8s2 â 1 â¤ 0, s1 + 2s2 + 11 â¥ 0}R4 = {s : s1 + 2s2 + 11 â¤ 0, s1 â¥ â3}R5 = {s : s1 â¤ â3, 3s1 + 4s2 + 25 â¤ 0}R6 = {s : 3s1 + 4s2 + 25 â¥ 0, 3s1 + 4s2 + 1 â¤ 0,â13s1 + 4s2 â 23 â¥ 0,C(s) â¥ 0}for C(s) = 9s21 + 24s1s2 â 234s1 + 16s22 + 200s2 â 527. It is shown inFigure 5.8.95Chapter 6Conclusions and Future workWe provided an algorithmic design to compute the closed convex enve-lope of a bivariate nonconvex quadratic function defined over a polytope.From [Loc16], we reformulated the original problem following the Assump-tion 3.2.1 to a subset of optimization subproblems. The solutions to thesesubproblems have either rational, quadratic or linear forms and are definedover R2 with a polyhedral subdivision. After solving the subproblems, theclosed convex envelope is computed as the maximum of all such solutions,which involves solving the map overlay problem over the original polytopefor all the polyhedral subdivisions of the solutionâs domains.The conjugate of the above rational function defined over a polytope Phas a parabolic subdivision such that over each member of the subdivisionit has form (4.13) - (4.15). The subdifferentials of the rational function isobserved to be contained in a parabola and when it is computed for int(P ) âR2, they lie along a parabolic arc. Similarly, for a quadratic function ofform (3.19) defined over a polytope, the conjugate is a piecewise definedfunction with either quadratic or linear form over a polyhedral subdivisionof R2.Unlike the computation of the closed convex envelope of a bilinear func-tion over a polytope, the conjugate has a hybrid data structure involvingboth the symbolic expressions and numerical values for the conjugate do-main, while the expressions are stored entirely in symbolic form. This iswell illustrated in Section 5.1.We formulated an algorithmic approach to compute the conjugate of therational functions of form (3.18). We start with computing the subdiffer-entials for the three entities of P , that is, vertices, edges and int(P ). Inthe special case of a vertex v where the function value is obtained by con-tinuity, the computation of subdifferentials is done at the end as âf(v) =R2 \u00E2xâPv2 âf(x). This region has been observed to be a union of polyhe-dral and parabolic region, and has been presented as Conjecture 4.11. If theconjecture does not hold, the conjugate would still be a fractional functiondefined over a parabolic subdivision of R2.966.1. Future workTable 6.1: Observed intersections for computing the maximum of all conju-gatesFunctions Linear Quadratic FractionalLinear line parabola lines or parabolaQuadratic parabola two lines lines or parabolaFractional lines or parabola lines or parabola lines or parabola6.1 Future workThe symbolic algorithm presented in Chapter 5.2 relies heavily on thesymbolic computations performed by the symbolic engine of the softwareMatlab. Compared to any numeric computation, symbolic computations arerelatively slow and a new module for computing the parabolic computationscould be created. However, similar to the polyhedral manipulations module,it will take significant developer time to do so.6.1.1 Maximum of all conjugatesThis thesis covers the computation of the conjugate of a single ratio-nal function defined over a polytope. However, the closed convex envelopeof a bivariate function is a piecewise rational function with a polyhedralsubdivision, so computing the conjugate of the piecewise function requirescomputing the conjugate of each piece, then finding the maximum of all suchconjugate expressions over each resulting parabolic subdivision. Now, sincewe are computing the maximum of the fractional functions of form (4.13),the points where two of them are equal may not lie on a parabolic arc (orparabola). However, for all the experiments conducted with valid conju-gates of a piecewise rational function which is the closed convex envelopeof a bilinear function, it has been observed to be parabolic. Table 6.1.1summarizes the observations.A parabolic intersection is also expected since taking the conjugate ofthe conjugate gives us the closure of our original convex piecewise ratio-nal function and by primal-dual symmetry, rational forms with polyhedraldomain lead to fractional form with parabolic domain, so the conjugate offractional forms will lead to rational form with polyhedral subdivision aswell.Conjecture 6.1. For a bivariate nonconvex quadratic function Q definedover a polytope P , let f(x) = conv QP (x) be a piecewise rational function976.1. Future workwith form (3.18)-(3.20), then the conjugate fâ(s) has a parabolic subdivi-sion such that over each member of its subdivision, it has one of the formsin (4.13)-(4.15).The computation of the maximum of all conjugates has been imple-mented in the module max conjugate, however, since we need to solve themap overlay problem for parabolic regions, which can increase exponentiallywith each added layer, it takes significant amount of computer time to returnthe resulting subdivision as the number of layers increases.6.1.2 Conjugate of bivariate nonconvex PLQ FunctionsIn the case of PLQ functions, we start with computing the closed convexenvelope of each piece and then compute the maximum of the conjugatesof each closed convex envelope piece to obtain the conjugate of the PLQfunction. Now, the conjugate is the maximum of all piecewise fractionalfunctions with parabolic subdivisions and since we do not know the struc-ture of the closed convex envelope of the PLQ function, we cannot inferanything about the resulting subdivision of the conjugate. From [Loc16],we have that the closed convex envelope of any bivariate function satisfyingAssumption 3.2.1 has a polyhedral subdivision, however, the form can onlybe implicitly defined. So if the closed convex envelope of a PLQ function hasrational form, then by symmetry we would know that the conjugate wouldhave a parabolic subdivision with fractional forms as well.Conjecture 6.2. For a bivariate nonconvex PLQ function g defined overa polytope P , let f(x) = conv gP (x), then if f is piecewise defined with apolyhedral subdivision of P such that over each member of its subdivisionit has form (3.18)-(3.20), then the conjugate fâ(s) has a parabolic subdivi-sion such that over each member of its subdivision, it has one of the formsin (4.13)-(4.15).6.1.3 Convex envelope of bivariate nonconvex PLQfunctionsMost of the PLQ functions defined over a polytope satisfy Assump-tion 3.2.1, so the closed convex envelope can be computed directly withthe method proposed in [Loc16]. For a quadratic function Q(x, y) = q1x2 +q2xy+ q3y2 + q4x+ q5y+ q6 and a convex edge of the polytope with relation986.1. Future worky = mjx+ qj , the η function is defined asâj â Eâ²(P ),ηj(a, b) =ï£±ï£²ï£³â(a+mjbâ Î³1j)24Î³2jâ bqj + γ0j (a, b) â Djâaxj(a, b)â byj(a, b) + xj(a, b)yj(a, b) (a, b) â {Dâj , D+j },where γ2j = q1 + q2mj + q3m2j , γ1j = q4 + q5mj + q2qj + 2mjq3qj andγ0j = q6 + q5qj + q3q2j . Now when solving the subproblems with quadratic-quadratic case, the relation between a and b is not linear, and so the onedimension optimization problems reduce to a quartic optimization problem.Similarly, the quadratic-linear case leads to a quartic optimization problemas well. So the solution may or may not have a rational form. However,if the solution has the same rational form as (3.18), the conjugate of PLQfunction would have a parabolic subdivision with the same fractional forms.Conjecture 6.3. For a bivariate nonconvex PLQ function g defined overa polytope P , if the solution of the one dimension optimization subprob-lems (3.17) has form (3.18)-(3.20), then conv gP has the same form and isdefined over a polyhedral subdivision of the original polytope P .99Bibliography[AB10] Kurt M Anstreicher and Samuel Burer. Computable representa-tions for convex hulls of low-dimensional quadratic forms. Math-ematical Programming, 124(1-2):33â43, 2010. â pages 4[AKF83] Faiz A Al-Khayyal and James E Falk. Jointly constrainedbiconvex programming. Mathematics of Operations Research,8(2):273â286, 1983. â pages 4[Ans12] Kurt M Anstreicher. On convex relaxations for quadraticallyconstrained quadratic programming. Mathematical Program-ming, 136(2):233â251, 2012. â pages 4[BGLW08] Heinz H Bauschke, Rafal Goebel, Yves Lucet, and Xianfu Wang.The proximal average: basic theory. SIAM Journal on Optimiza-tion, 19(2):766â785, 2008. â pages 1[BLT08] Heinz H Bauschke, Yves Lucet, and Michael Trienis. How totransform one convex function continuously into another. SIAMReview, 50(1):115â132, 2008. â pages 1, 2[Bre89] Yann Brenier. Un algorithme rapide pour le calcul de trans-forme´es de Legendre-Fenchel discretes. Comptes rendus delâAcade´mie des sciences. Se´rie 1, Mathe´matique, 308(20):587â589, 1989. â pages 1, 4[CCHP17] Ying Cui, Tsung-Hui Chang, Mingyi Hong, and Jong-Shi Pang.A study of piecewise linear-quadratic programs. arXiv preprintarXiv:1709.05758, 2017. â pages 1[Cor96] Lucilla Corrias. Fast Legendre-Fenchel transform and applica-tions to Hamilton-Jacobi equations and conservation laws. SIAMJournal on Numerical Analysis, 33(4):1534â1558, 1996. â pages1, 4100Bibliography[Cra89] Yves Crama. Recognition problems for special classes of polyno-mials in 0â1 variables. Mathematical Programming, 44(1-3):139â155, 1989. â pages 3[GJL14] Bryan Gardiner, Khan Jakee, and Yves Lucet. Computing thepartial conjugate of convex piecewise linear-quadratic bivari-ate functions. Computational Optimization and Applications,58(1):249â272, 2014. â pages 5[GL10] Bryan Gardiner and Yves Lucet. Convex hull algorithms forpiecewise linear-quadratic functions in computational convexanalysis. Set-Valued and Variational Analysis, 18(3-4):467â482,2010. â pages 2[GL11] Bryan Gardiner and Yves Lucet. Graph-matrix calculus forcomputational convex analysis. In Fixed-Point Algorithms forInverse Problems in Science and Engineering, pages 243â259.Springer, 2011. â pages 1[GL13] Bryan Gardiner and Yves Lucet. Computing the conjugate ofconvex piecewise linear-quadratic bivariate functions. Mathe-matical Programming, 139(1-2):161â184, 2013. â pages 5, 59[Her16] Cristopher Hermosilla. Legendre transform and applicationsto finite and infinite optimization. Set-Valued and VariationalAnalysis, 24(4):685â705, 2016. â pages 4[HL18] Tasnuva Haque and Yves Lucet. A linear-time algorithm to com-pute the conjugate of convex piecewise linear-quadratic bivari-ate functions. Computational Optimization and Applications,70(2):593â613, 2018. â pages 5[HUL93] Jean-Baptiste Hiriart-Urruty and Claude Lemare´chal. Convexanalysis and minimization algorithms II: Advanced Theory andBundle Methods. Springer Science & Business Media, 1993. âpages 5, 6[HUL07] Jean-Baptiste Hiriart-Urruty and Yves Lucet. Parametric com-putation of the Legendre-Fenchel conjugate with application tothe computation of the Moreau envelope. Journal of ConvexAnalysis, 14(3):657, 2007. â pages 1101Bibliography[HUL13] Jean-Baptiste Hiriart-Urruty and Claude Lemare´chal. Convexanalysis and minimization algorithms I: Fundamentals, volume305. Springer Science & Business Media, 2013. â pages 8, 11,14, 15, 50[JKL11] Jennifer A Johnstone, Valentin R Koch, and Yves Lucet. Con-vexity of the proximal average. Journal of Optimization Theoryand Applications, 148(1):107â124, 2011. â pages 1[JMW08] Matthias Jach, Dennis Michaels, and Robert Weismantel. Theconvex envelope of (nâ1)-convex functions. SIAM Journal onOptimization, 19(3):1451â1466, 2008. â pages 3[KBB+09] Leonid Khachiyan, Endre Boros, Konrad Borys, Vladimir Gur-vich, and Khaled Elbassioni. Generating all vertices of a poly-hedron is hard. In Twentieth Anniversary Volume:, pages 1â17.Springer, 2009. â pages 48[LBT09] Yves Lucet, Heinz H Bauschke, and Mike Trienis. The piecewiselinear-quadratic model for computational convex analysis. Com-putational Optimization and Applications, 43(1):95â118, 2009.â pages 1, 2[Lin86] Nathan Linial. Hard enumeration problems in geometry andcombinatorics. SIAM Journal on Algebraic Discrete Methods,7(2):331â335, 1986. â pages 48[Lin05] Jeff Linderoth. A simplicial branch-and-bound algorithm forsolving quadratically constrained quadratic programs. Mathe-matical Programming, 103(2):251â282, 2005. â pages 4[LLT17] Florian Lauster, D Russell Luke, and Matthew K Tam. Symboliccomputation with monotone operators. Set-Valued and Varia-tional Analysis, pages 1â16, 2017. â pages 1[Loc14] Marco Locatelli. A technique to derive the analytical form ofconvex envelopes for some bivariate functions. Journal of GlobalOptimization, 59(2-3):477â501, 2014. â pages 4[Loc16] Marco Locatelli. Polyhedral subdivisions and functional formsfor the convex envelopes of bilinear, fractional and other bivari-ate functions over general polytopes. Journal of Global Opti-mization, 66(4):629â668, 2016. â pages v, vi, 4, 7, 13, 14, 16,17, 18, 19, 20, 21, 25, 26, 27, 28, 29, 96, 98102Bibliography[LS14] Marco Locatelli and Fabio Schoen. On convex envelopes forbivariate functions over polytopes. Mathematical Programming,144(1-2):65â91, 2014. â pages 18[Luc96] Yves Lucet. A fast computational algorithm for the Legendre-Fenchel transform. Computational Optimization and Applica-tions, 6(1):27â57, 1996. â pages 1, 4[Luc97] Yves Lucet. Faster than the fast Legendre transform, the linear-time Legendre transform. Numerical Algorithms, 16(2):171â185,1997. â pages 1, 4[Luc06] Yves Lucet. Fast Moreau envelope computation i: Numerical al-gorithms. Numerical Algorithms, 43(3):235â249, 2006. â pages1[Luc10] Yves Lucet. What shape is your conjugate? a survey of com-putational convex analysis and its applications. SIAM Review,52(3):505â542, 2010. â pages 1, 2[McC76] Garth P McCormick. Computability of global solutions to fac-torable nonconvex programs: Part iâconvex underestimatingproblems. Mathematical Programming, 10(1):147â175, 1976. âpages 4[Mor65] Jean-Jacques Moreau. Proximite´ et dualite´ dans un espacehilbertien. Bull. Soc. Math. France, 93(2):273â299, 1965. âpages 1[NV94] A Noullez and M Vergassola. A fast Legendre transform al-gorithm and applications to the adhesion model. Journal ofScientific Computing, 9(3):259â281, 1994. â pages 1[RW98] R Tyrrell Rockafellar and Roger J-B Wets. Variational Analysis,volume 317. Springer Science & Business Media, 1998. â pages1[SA90] Hanif D Sherali and Amine Alameddine. An explicit characteri-zation of the convex envelope of a bivariate bilinear function overspecial polytopes. Annals of Operations Research, 25(1):197â209, 1990. â pages 4103Bibliography[SAF92] Zhen-Su She, Erik Aurell, and Uriel Frisch. The inviscid Burgersequation with initial data of brownian type. Communications inMathematical Physics, 148(3):623â641, 1992. â pages 1[Tar04] Fabio Tardella. On the existence of polyhedral convex envelopes.In Frontiers in Global Optimization, pages 563â573. Springer,2004. â pages 3[Tar08] Fabio Tardella. Existence and sum decomposition of vertex poly-hedral convex envelopes. Optimization Letters, 2(3):363â375,2008. â pages 3104Appendix105Appendix AThe η(a, b) expressions for thebilinear functionsLet f(x, y) = xy be defined over a polytope P . Now the restriction off along the edge between vertices vj1 and vj2 , defined by points statisfyingyâ² = mjxâ² + qj , isfj(xâ²) = xâ² · (mjxâ² + qj)= mjxâ²2 + qjxâ²,and its Hessianâ2fj(xâ²) = 2mj .So fj can only be strictly convex when mj > 0.Now the minimization problem related to this edge can be written asminxâ²â[xvj1 ,xvj2 ]fj(xâ²)â [a(xâ² â x) + b(yâ² â y)] â¥ cmjxâ²2 + qjxâ² â (a+ bmj)xâ² â bqj + ax+ by â¥ c (A.1)Let C(xâ²) = mjxâ²2 + qjxâ² â (a + bmj)xâ² + ax + by â bqj . Then theoptimal value of the above unconstrained problem, obtained by computingthe critical points, is2mjxâ² + qj â (a+ bmj) = 0xâ² =a+ bmj â qj2mj.So the optimal value of (A.1) isxj(a, b) =ï£±ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£³xvj1 (a, b) â Dâjxvj2 (a, b) â D+ja+ bmj â qj2mj(a, b) â Dj106Appendix A. The η(a, b) expressions for the bilinear functionswhere the three regions (defined in (3.11)) areDâj = {(a, b) : a+ bmj â¤ 2mjxvj1 + qj}D+j = {(a, b) : a+ bmj â¥ 2mjxvj2 + qj}Dj = {(a, b) : 2mjxvj1 + qj â¤ a+ bmj â¤ 2mjxvj2 + qj}.Also, yj(a, b) = mjxj(a, b) + qj . Now from (3.13), âk â V â²(P ) âª Eâ²(P )(see(3.8)),ηk(a, b) = f(xk(a, b), yk(a, b))â axk(a, b)â byk(a, b).So for f(x, y) = xy and j â Eâ²(P ),ηj(a, b) = xj(a, b)yj(a, b)â axj(a, b)â byj(a, b).Now when (a, b) â Dj ,ηj(a, b) = mj(a+ bmj â qj2mj)2+ qja+ bmj â qj2mjâ aa+ bmj â qj2mjâ b(mja+ bmj â qj2mj+ qj)=(a+ bmj â qj2mj)·(a+ bmj â qj2+ qj â aâ bmj)â bqj=(a+ bmj â qj4mj)· (a+ bmj â qj + 2qj â 2aâ 2bmj)â bqj=(a+ bmj â qj4mj)· (âaâ bmj + qj)â bqj=â (a+ bmj â qj)24mjâ bqj .So âj â Eâ²(P ),ηj(a, b) =ï£±ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£³âaxvj1 â byvj1 + xvj1yvj1 (a, b) â Dâjâ(a+mjbâ qj)24mjâ bqj (a, b) â Djâaxvj2 â byvj2 + xvj2yvj2 (a, b) â D+j .Similarly, âi â V â²(P ),ηi(a, b) = xviyvi â axvi â byvi .107Appendix BSolving the subproblemsFor h,w â Eâ²(P ), over a subregion Sr, the subproblem ismaxa,bηh(a, b) + ax+ byηh(a, b) = ηw(a, b)ηh(a, b) â¤ ηk(a, b), k â Eâ²(P ) âª V â²(P ) {h,w}(a, b) â SrSince for the bilinear term, η is either quadratic or linear, solving thesubproblems lead to the three cases: quadratic-quadratic, quadratic-linearand linear-linear as discussed in Section 3.2.3.B.1 Case 1: Quadratic-QuadraticLet the respective convex edges for ηh and ηw be characterized by mh >0, qh and mw > 0, qw. Soηh(a, b) = â(a+ bmh â qh)24mhâ bqh, and ηw(a, b) = â(a+ bmw â qw)24mwâ bqw.Now, in this case, the subproblem ismaxa,bâ (a+ bmh â qh)24mhâ bqh + ax+ byâ (a+ bmh â qh)24mhâ bqh = â(a+ bmw â qw)24mwâ bqwâ (a+ bmh â qh)24mhâ bqh â¤ ηk(a, b), k â Eâ²(P ) âª V â²(P ) {h,w}lh â¤ a+ bmh â¤ uh, lw â¤ a+ bmw â¤ uw, · · Â·where lh â¤ a+ bmh â¤ uh, lw â¤ a+ bmw â¤ uw, · · · are the active constraintsfor the subregion Sr.108B.1. Case 1: Quadratic-QuadraticB.1.1 Solving the equality constraintNow we start with solving the equality constraint asâ(a+ bmh â qh)24mhâ bqh = â(a+ bmw â qw)24mwâ bqwBy rearranging the terms, we obtain(mh âmw)a2 â 2(mhqw âmwqh)a+mh(mwb+ qw)2 âmw(mhb+ qh)2 = 0.When mh = mw,â 2mh(qw â qh)a+mh((mhb)2 + 2mhqwb+ q2w)âmh((mhb)2 + 2mhqhb+ q2h) = 0,or â 2mh(qw â qh)a+ 2m2hqwbâ 2m2hqhb+mh(q2w â q2h) = 0,or â 2(qw â qh)a+ 2mh(qw â qh) + q2w â q2h = 0,or â 2(qw â qh)a+ 2mh(qw â qh) + (qw + qh) · (qw â qh) = 0,or â 2a+ 2mhb+ (qh + qw) = 0,or a = mhb+qh + qw2,and when mh 6= mw,a =mhqw â qhmw Â±âmhmw[(mh âmw)b+ qh â qw]mh âmw .Now using the above relation(s), we reduce the subregion {(a, b) : lh â¤a+ bmh â¤ uh, lw â¤ a+ bmw â¤ uw, · · Â·} to simple bound constraints asl â¤ b â¤ u.If the intersection of all such one dimension bounds for all the active con-straints is an empty set, the subproblem becomes infeasible and we proceedon to solving the next pair.B.1.2 Solving the inequalitiesWhen the one dimesion subregion is feasible, using the relations, we solvethe inequalities for all the remaining Î·k,â k â Eâ²(P )âªV â²(P ){h,w}, whichleads to two cases: quadratic-quadratic and quadratic-linear. For simplicity,let the relation be denoted by a = α1b+ α0.109B.1. Case 1: Quadratic-QuadraticCase 1.1 : Quadratic-Quadratic inequalityWhen ηk is quadratic, the inequality constraint isâ (a+ bmh â qh)24mhâ bqh â¤ â(a+ bmk â qk)24mkâ bqkor â ((α1 +mh)b+ α0 â qh)24mhâ bqh â¤ â((Î±1 +mk)b+ α0 â qk)24mkâ bqkUpon simplification, it reduces to the formβ2b2 + β1b+ β0 â¤ 0whereβ2 = mh((α1 +mk)2)âmk((Î±1 +mh)2),β1 = 2mh(α1 +mk) · (α0 â qk)â 2mk(α1 +mh) · (α0 â qh) + 4mhmk(qkâqh),Î²0 = mh((α0 â qk)2)âmk((Î±0 â qh)2).When β2 = 0, the inequality reduces to an interval, otherwise, the solutionis a union of intervals.Case 1.2 : Quadratic-Linear inequalityWhen ηk is linear, the inequality constaint isâ (a+ bmh â qh)24mhâ bqh â¤ âxkaâ ykb+ xkyk,or â ((α1 +mh)b+ α0 â qh)24mhâ bqh â¤ âxk(Î±1b+ Î±0)â ykb+ xkyk.Upon simplication, it reduces toβ2b2 + β1b+ β0 â¤ 0withβ2 = â(Î±1 +mh)24mhβ1 = â(Î±1 +mh) · (α0 â qh)2mh+ yk + xkα1 â qhβ0 = â(Î±0 â qh)24mh+ xkα0 â xkyk.110B.1. Case 1: Quadratic-QuadraticSimilar to Case 1.1 in Section B.1.2, when β2 = 0, the inequality reduces toan interval, otherwise, to a union of intervals.We solve the inequalities for all ηk, at the same time, verify that thesolution, along with the intersection with the one dimension subregion, isfeasible.After solving all the inequalities, the bounds on b reduce to the formdl â¤ b â¤ du.B.1.3 SolutionsNext, using the relation, we reduce the subproblem to one dimension asmax â ((α1 +mh)b+ α0 â qh)24mhâ bqh + (α1b+ α0)x+ bydl â¤ b â¤ du.By simplifying the terms, it reduces tomax â ξ2b2 + 2bξ1(x, y) + ξ0(x, y) (B.1)dl â¤ b â¤ duwhereξ2 =(α1 +mh)24mhâ¥ 0,ξ1(x, y) = â(Î±1 +mh)(α0 â qh)2mh+ y â qh + α1x,ξ0(x, y) = â(Î±0 â qh)24mh+ α0xas mh > 0.Since (B.1) is concave and differentiable in b, the optimal solution of theabove unconstrained problem isddb(âÎ¾2b2 + 2bξ1(x, y) + ξ0(x, y)) = 0,or â 2ξ2b+ 2ξ1(x, y) = 0,or b =ξ1(x, y)ξ2.111B.2. Case 2: Quadratic-LinearConsequently, its optimal value= âÎ¾2(Î¾1(x, y)ξ2)2 + 2ξ1(x, y)ξ2ξ1(x, y) + ξ0(x, y)= âÎ¾21(x, y)ξ2+ 2ξ21(x, y)ξ2+ ξ0(x, y)=ξ21(x, y)ξ2+ ξ0(x, y).So the optimal solution of the subproblem (B.1) isï£±ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£³dl whenÎ¾1(x,y)Î¾2â¤ dlξ1(x, y)ξ2when dl <ξ1(x,y)ξ2< dudu otherwise,and hence its optimal value isï£±ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£³âÎ¾2d2l + 2dlξ1(x, y) + ξ0(x, y) if dlξ2 â¥ ξ1(x, y),ξ21(x, y)ξ2+ ξ0(x, y) if dlξ2 < ξ1(x, y) < duÎ¾2,âÎ¾2d2u + 2duξ1(x, y) + ξ0(x, y) otherwise.B.2 Case 2: Quadratic-LinearLet the respective convex edge for ηh be characterized by mh > 0, qh andthe vertex for ηw by (xw, yw), thenηh = â(a+ bmh â qh)24mhâ bqh and ηw = âxwaâ ywb+ xwyw.So the subproblem ismaxa,bâ (a+ bmh â qh)24mhâ bqh + ax+ byâ (a+ bmh â qh)24mhâ bqh = âxwaâ ywb+ xwywâ (a+ bmh â qh)24mhâ bqh â¤ ηk(a, b), k â Eâ²(P ) âª V â²(P ) {h,w}lh â¤ a+ bmh â¤ uh, · · Â·112B.2. Case 2: Quadratic-LinearB.2.1 Solving the equality constraintSince the equality constraint isâ(a+ bmh â qh)24mhâ bqh = âxwaâ ywb+ xwyw,let zh = a+mhb, then a = zh âmhb. Soâ(zh â qh)24mhâ bqh = âxw(zh âmhb)â ywb+ xwyw,which upon rearranging results inb =z2h â (4mhxw + 2qh)zh + 4mhxwyw + q2h4mh(mhxw + q â yw) .For simplicity, let b = β2z2h + β1zh + β0. So the equality constraint resultsin a quadratic relation for a and b in zh.Now using the above relations, we reduce the subregion {(a, b) : lh â¤a+ bmh â¤ uh, · · Â·} to one dimension intervals asl â¤ zh â¤ u.Similar to Case 1 quadratic-quadratic , when the intersection results in anempty set, the subproblem is declared infeasible and we proceed on to solvingthe next pair.B.2.2 Solving the inequalitiesSolving the inequalities for all the remaining ηk defined on the subregion,similar to Case 1(see SectionB.1.2), leads to two cases as quadratic-quadraticand quadratic-linear.Case 2.1 : Quadratic-Quadratic inequalityWhen ηk is quadratic, the inequality constraint isâ (zh â qh)24mhâ (β2z2h + β1zh + β0)qh â¤ â(a+ bmk â qk)24mkâ bqkorâ (zh â qh)24mhâ (β2z2h + β1zh + β0)qh â¤â (zh + (mk âmh)(Î²2z2h + β1zh + Î²0)â qk)24mkâ (β2z2h + β1zh + β0)qk.113B.2. Case 2: Quadratic-LinearUpon simplication, it reduces to a quartic inequality asζ4z4h + ζ3z3h + ζ2z2h + ζ1zh + ζ0 â¤ 0whereζ4 =t224mkζ3 =2t2t14mkζ2 =t21 + 2t2t04mk+ β2(qk â qh)â 14mhζ1 =2t1t04mk+ β1(qk â qh) + 2qh4mhζ0 =t204mk+ β0(qk â qh)â q2h4mh,and t2 = Î²2(mkâmh), t1 = Î²1(mkâmh)+1, t0 = Î²0(mkâmh)âqk. Whenmh = mk, the inequality reduces to quadratic form.The solution, in this case, is a union of intervals for zh.Case 2.2 : Quadratic-Linear inequalityWhen ηk is linear, the inequality constraint isâ (zh â qh)24mhâ (β2z2h + β1zh + β0)qh â¤â xkzh + (xkmh â yk) · (β2z2h + β1zh + β0) + xkyk.By rearranging the terms, we getζ2z2h + ζ1zh + ζ0 â¤ 0whereζ2 = β2(yk +mhxk â qh)â 14mhζ1 = β1(yk + xkmh â qh)â xk + qh2mhζ0 = β0(yk +mhxk â qh)â xkyk â q2h4mh.The solution of the above quadratic inequality is again a union of the inter-vals for zh.114B.2. Case 2: Quadratic-LinearB.2.3 SolutionsBy subtituting the expressions for a and b, the subproblem becomesmax â (a+ bmh â qh)24mhâ bqh + ax+ bydl â¤ a+ bmh â¤ du= max â (zh â qh)24mhâ (β2z2h + β1zh + β0)qh + zhx+ (y â xmh)·(β2z2h + β1zh + β0)dl â¤ zh â¤ du.By rearranging the terms, it reduces to the formmax â ξ2(x, y)z2h + 2zhξ1(x, y) + ξ0(x, y) (B.2)dl â¤ zh â¤ duwhereξ2(x, y) = Î²2xâ β2y + β2qh + 14mh,ξ1(x, y) =12[(1â β1)x+ β1y â β1qh + qh2mh],ξ0(x, y) = âÎ²0x+ β0y â β0qh â q2h4mh.Now when ξ2(x, y) > 0, the optimal solution for the unconstrained prob-lem of (B.2) isddzh(âÎ¾2(x, y)z2h + 2zhξ1(x, y) + ξ0(x, y)) = 0or â 2ξ2(x, y)zh + 2ξ1(x, y) = 0or zh =ξ1(x, y)ξ2(x, y)and consequently, its optimal value= âÎ¾2(x, y)(ξ1(x, y)ξ2(x, y))2 + 2ξ1(x, y)ξ2(x, y)ξ1(x, y) + ξ0(x, y)= âÎ¾21(x, y)ξ2(x, y)+ 2ξ21(x, y)ξ2(x, y)+ ξ0(x, y)=ξ21(x, y)ξ2(x, y)+ ξ0(x, y).115B.3. Case 3: Linear-LinearSo the optimal solution of (B.2), when ξ2(x, y) > 0 isï£±ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£³dl whenÎ¾1(x,y)Î¾2(x,y)â¤ dlξ1(x, y)ξ2(x, y)when dl <ξ1(x,y)ξ2(x,y)< dudu otherwiseand its optimal value isï£±ï£´ï£´ï£´ï£²ï£´ï£´ï£´ï£³âÎ¾2(x, y)d2l + 2dlξ1(x, y) + ξ0(x, y) if dlξ2(x, y) â¥ ξ1(x, y),ξ21(x, y)ξ2(x, y)+ ξ0(x, y) if dlξ2(x, y) < ξ1(x, y) < duξ2(x, y),âÎ¾2(x, y)d2u + 2duξ1(x, y) + ξ0(x, y) otherwise,and when ξ2(x, y) â¤ 0, it ismax{ â ξ2(x, y)d2l + 2dlξ1(x, y) + ξ0(x, y),âÎ¾2(x, y)d2u + 2duξ1(x, y)+ξ0(x, y)}.B.3 Case 3: Linear-LinearLet the two vertices corresponding to ηh and ηw be defined as (xh, yh)and (xw, yw). Then the subproblem in this case is defined asmaxa,bâ xhaâ yhb+ xhyh + ax+ byâ xhaâ yhb+ xhyh = âxwaâ ywb+ xwywâ xhaâ yhb+ xhyh â¤ ηk(a, b), k â Eâ²(P ) âª V â²(P ) {h,w}lk â¤ a+ bmk â¤ uk, · · Â·B.3.1 Solving the equality constraintWe start with solving the equality constraint first. Since both ηh and ηware linear, the relationâ xhaâ yhb+ xhyh = âxwaâ ywb+ xwywor a =yh â ywxw â xh b+xhyh â xwywxh â xwis linear as well. For simplicity, let a = α1b+ α0.By using this relation, we reduce the subregion to simple one-dimensionbounds as l â¤ b â¤ u. When the intersection results in an empty set, thesubproblem is infeasible.116B.3. Case 3: Linear-LinearB.3.2 Solving the inequalitiesNext we solve the inequalities for the remaining ηk.Case 3.1 : Linear-QuadraticLet the convex edge corresponding to ηk be characterized by mk and qk.Thenâ xhaâ yhb+ xhyh â¤ â(a+ bmk â qk)24mkâ bqk,or â xh(α1b+ Î±0)â yhb+ xhyh â¤ â((Î±1 +mk)b+ α0 â qk)24mkâ bqk.Upon simplification this leads to a quadratic inequality in b asβ2b2 + β1b+ β0 â¤ 0withβ2 =(α1 +mk)24mkβ1 =(α1 +mk) · (α0 â qk)2mk+ qk â α1xh â yhβ0 =(α0 â qk)24mk+ xhyh â xhα0.So the solution in this case is a union of intervals.Case 3.2 : Linear-LinearWhen ηk is linear, the inequality reduces toâ xhaâ yhb+ xhyh â¤ âxkaâ ykb+ xkyk,or â xh(α1b+ Î±0)â yhb+ xhyh â¤ âxk(Î±1b+ Î±0)â ykb+ xkyk,or (α1xk + yk â α1xh â yh)b â¤ xkyk + α0xh â xhyh â α0xk.Now when α1xk + yk â α1xh â yh > 0,b â¤ xkyk + α0xh â xhyh â α0xkα1xk + yk â α1xh â yhand when α1xk + yk â α1xh â yh < 0,b â¥ xkyk + α0xh â xhyh â α0xkα1xk + yk â α1xh â yh .117B.3. Case 3: Linear-LinearB.3.3 SolutionsBy substituting the linear relation between a and b in the subproblem,we getmaxa,bâ xh(α1b+ Î±0)â yhb+ xhyh + (α1 + α0)x+ bydl â¤ b â¤ du.Upon simplification, it reduces tomax ξ1(x, y)b+ ξ0(x, y)dl â¤ b â¤ duwhereξ1(x, y) = α1x+ y â α1xh â yh,ξ0(x, y) = α0x+ xhyh â α0xh.Since we are maximizing a linear function, the optimal value lies at theboundaries, that is,max {dlξ1(x, y) + ξ0(x, y), duξ1(x, y) + ξ0(x, y)}.118Appendix CConjugate expressions for arational function over anedgeFor a bivariate rational function r(x) =ξ21(x)ξ2(x)+ ξ0(x) with ξ2(x) =ξ21x1+ξ22x2+ξ20, ξ1(x) = ξ11x1+ξ12x2+ξ10 and ξ0(x) = ξ01x1+ξ02x2+ξ00,and a polytope P , let f(x) = r(x)+IP (x) and E = {x : x2 = mx1+q, vâ1 â¤x1 â¤ v+1} â P be an edge between vertices vâ and v+. Then the conjugatefâ(s) = supxâri(E){ãs, xã â f(x)}= supxâri(E){ãs, xã â (r(x) + IP (x))}.By computing the critical points, we getsâ (âr(x) +NP (x)) = 0.Now NP (x) = {s : s = Î»(âm, 1), λ â¥ 0}, sos1 = âÎ¾21t2 + 2ξ11t+ ξ01 âmÎ»s2 = âÎ¾22t2 + 2ξ12t+ ξ02 + λ (C.1)where t =ξ11x1 + ξ12x2 + ξ10ξ21x1 + ξ22x2 + ξ20. As x â ri(E), we also havex2 = mx1 + q. (C.2)Now by eliminating λ from (C.1), we obtains1 +ms2 = â(Î¾21 +mξ22)t2 + 2(ξ11 +mξ12)t+ ξ01 +mξ02). (C.3)119Appendix C. Conjugate expressions for a rational function over an edgeWhen ξ21 +mξ22 = 0, (C.3) becomess1 +ms2 = 2(ξ11 +mξ12)t+ ξ01 +mξ02,or t =s1 +ms2 â ξ01 âmÎ¾022(Î¾11 +mξ12).Let t = t1s1 + t2s2 + t0, then by using the value of t and (C.2), we getξ11x1 + ξ12x2 + ξ10ξ21x1 + ξ22x2 + ξ20= t1s1 + t2s2 + t0or ξ11x1 + (q +mx1)ξ12 + ξ10 â (t1s1 + t2s2 + t0) · (ξ20 + ξ21x1 + (q+mx1)ξ22) = 0or x1 =(ξ20 + qξ22) · (t0 + t1s1 + t2s2)â ξ10 â qξ12ξ11 +mξ12.For simplification, let x1 = γ10s1 + γ01s2 + γ00 whereγ10 =t1(ξ20 + qξ22)ξ11 +mξ12γ01 =t2(ξ20 + qξ22)ξ11 +mξ12γ00 =t0(ξ20 + qÎ¾22)â ξ10 â qξ12ξ11 +mξ12,then the conjugate expression, upon substituting values, becomesfâ(s) = s1x1 + s2(mx1 + q)â r(x1,mx1 + q)= ζ11s21 + ζ12s1s2 + ζ22s22 + ζ10s1 + ζ01s2 + ζ00whereζ11 = â(Î¾11Î³10 +mξ12γ10)2ξ20 + qξ22+ γ10ζ12 = â2(Î¾11Î³01 +mξ12γ01) · (ξ11γ10 +mξ12γ10)ξ20 + qξ22+ γ01 +mγ10ζ22 = â(Î¾11Î³01 +mξ12γ01)2ξ20 + qξ22+ γ01mζ10 = â2(Î¾11Î³10 +mξ12γ10) · (ξ10 + ξ11γ00 + ξ12(q +mγ00))ξ20 + qÎ¾22âmÎ¾02Î³10+ γ00 â ξ01γ10120Appendix C. Conjugate expressions for a rational function over an edgeζ01 = â2(Î¾11Î³01 +mξ12γ01) · (ξ10 + ξ11γ00 + ξ12(q +mγ00))ξ20 + qÎ¾22âmÎ¾02Î³01â ξ01γ01 +mγ00 + qζ00 = â(Î¾10 + ξ11γ00 + ξ12(q +mγ00))2ξ20 + qÎ¾22â ξ00 â ξ01γ00 â ξ02(q +mγ00).Now when ξ21 +mξ22 6= 0, (C.3) results ins1 +ms2 = â(Î¾21 +mξ22)t2 + 2(ξ11 +mξ12)t+ ξ01 +mξ02,or (ξ21 +mξ22)t2 â 2(ξ11 +mξ12)t+ s1 +ms2 â ξ01 âmÎ¾02 = 0,or t =t00 Â±ât10/2s1 + t01/2s2 + t00/2tâ00wheret00 = ξ11 +mξ12t10/2 = âÎ¾21 âmÎ¾22t01/2 = âÎ¾22m2 âmÎ¾21t00/2 = ξ01ξ21 + ξ211 +m(mξ212 + ξ01ξ22 + ξ02ξ21 + 2ξ11ξ12 +mÎ¾02Î¾22)tâ00 = ξ21 +mξ22.For simplicity, let t1/2 =ât10/2s1 + t01/2s2 + t00/2. By using t =ξ11x1 + ξ12x2 + ξ10ξ21x1 + ξ22x2 + ξ20and (C.2), x2 = mx1 + q, we getξ11x1 + ξ12x2 + ξ10ξ21x1 + ξ22x2 + ξ20=t00 Â±ât10/2s1 + t01/2s2 + t00/2tâ00,or ξ10 + ξ11x1 + ξ12x2 â(t00 Â± t1/2) · (ξ20 + ξ21x1 + Î¾22x2)tâ00= 0,or x1 =ξ20t00 â Î¾10tâ00 Â± ξ20t1/2 â qÎ¾12tâ00 + qξ22t00 Â± qÎ¾22t1/2Î¾11tâ00 â ξ21t00 â ξ21t1/2 +mÎ¾12tâ00 âmÎ¾22t00 âmÎ¾22t1/2,or x1 =γ00 Â± Î³1/2âÎ³10/2s1 + γ01/2s2 + Î³00/2âÎ³â1/2âÎ³10/2s1 + γ01/2s2 + γ00/2121Appendix C. Conjugate expressions for a rational function over an edgewhereγ00 = ξ20t00 â Î¾10tâ00 â qÎ¾12tâ00 + qξ22t00γ10/2 = t10/2γ01/2 = t01/2γ00/2 = t00/2Î³â1/2 = âÎ¾21 âmÎ¾22.Now by substituting the optimal values of x1 and x2, the conjugateexpression becomesfâ(s) = s1x1 + s2(mx1 + q)â r(x1,mx1 + q)=Ï1(s1, s2)Î¶00âÏ1/2(s1, s2)+ Ï0(s1, s2)whereζ00 = (ξ21 +mÎ¾22)2Ï1(s1, s2) = 2(ξ21 +mξ22) · (ξ10ξ21 â ξ11ξ20 +mξ10ξ22 âmÎ¾12Î¾20 â qξ11ξ22+ qξ12ξ21)s1 + 2m(mξ22 + ξ21) · (ξ10ξ21 â ξ11ξ20 +mÎ¾10Î¾22âmÎ¾12Î¾20 â qξ11ξ22 + qξ12ξ21)s2 + Î´1Ï1/2(s1, s2) = â(Î¾21 +mξ22)s1 âm(Î¾21 +mξ22m)s2 + ξ01ξ21 + ξ211+m(mξ212 + ξ01ξ22 + ξ02ξ21 + 2ξ11ξ12 +mÎ¾02Î¾22)Ï0(s1, s2) = (âÎ¾20(Î¾21 +mÎ¾22)â qξ22(ξ21 +mξ22))s1 + (qξ21(ξ21 +mÎ¾22)âmÎ¾20(Î¾21 +mξ22))s2 + δ0.andδ1 =â 2(ξ10ξ21 â ξ11ξ20 +mξ10ξ22 âmÎ¾12Î¾20 â qξ11ξ22 + qξ12ξ21)· (ξ01ξ21 + ξ211 +m(ξ212m+ ξ01ξ22 + ξ02ξ21 + 2ξ11ξ12 + ξ02ξ22m)).δ0 = 2ξ11ξ20(ξ11 +mÎ¾12)â 2ξ10ξ21(ξ11 +mÎ¾12)â ξ00ξ21(ξ21 +mξ22)+ ξ01ξ20(ξ21 +mÎ¾22)â 2mξ10ξ22(ξ11 +mξ12) + 2mξ12ξ20(ξ11 +mÎ¾12)âmÎ¾00Î¾22(Î¾21 +mξ22) +mξ02ξ20(ξ21 +mξ22) + 2qξ11ξ22(ξ11 +mÎ¾12)â 2qξ12ξ21(ξ11 +mξ12) + qξ01ξ22(ξ21 +mÎ¾22)â qξ02ξ21(ξ21 +mξ22).Hence when ξ21 +mξ22 6= 0,fâ(s) =Ï1(s1, s2)Î¶00âÏ1/2(s1, s2)+ Ï0(s1, s2),122Appendix C. Conjugate expressions for a rational function over an edgeand when ξ21 +mξ22 = 0,fâ(s) = ζ11s21 + ζ12s1s2 + ζ22s22 + ζ10s1 + ζ01s2 + ζ00,where all ζij , Ïi and Ïi/j have been formulated above.123
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Conjugate of some rational functions and convex envelope of quadratic functions over a polytope
Text