Science, Irving K. Barber Faculty of (Okanagan)
Computer Science, Mathematics, Physics and Statistics, Department of (Okanagan)
DSpace
UBCO
Karmarkar, Tanmaya Narendra
Master of Science - MSc
University of British Columbia
Computing the minima of a convex function is much easier than computing the same for a nonconvex function. Every nonconvex function has the same minimum as its convex envelope and the minima of the nonconvex function are contained in the set of minima of its convex envelope. We compute the conjugate of a bivariate piecewise linear-quadratic (PLQ) function as a first step toward the computation of the convex envelope. Our algorithm starts with computing the convex envelope of each piece obtaining a rational function defined over a polyhedral subdivision. Then we compute the conjugate of each of those pieces and obtain a fractional form defined over a parabolic subdivision. The last step is to compute the maximum of all those functions to obtain the conjugate of the original piecewise linear-quadratic function as a piecewise function defined on a parabolic subdivision. Our implementation in MATLAB uses symbolic computation and rational numbers to avoid any floating-point errors. We have developed the theory to compute the conjugate by computing the maximum of the conjugates corresponding to all the pieces and have implemented the code to compute the convex envelope, conjugate of each piece and the conjugate of the PLQ function.
Computing the conjugate ofnonconvex bivariate piecewiselinear-quadratic functionsbyTanmaya Narendra KarmarkarB.E. Electrical Engineering, College of Engineering, Pune 1999A 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)July 2024© Tanmaya Narendra Karmarkar, 2024The following individuals certify that they have read, and recommendto the College of Graduate Studies for acceptance, a thesis/dissertation en-titled:Computing the conjugate of nonconvex bivariate piecewiselinear-quadratic functionssubmitted by Tanmaya Narendra Karmarkar in partial fulfilment ofthe requirements of the degree of Master of ScienceYves Lucet, I. K. Barber School Faculty of ScienceSupervisorHeinz Bauschke, I. K. Barber School Faculty of ScienceSupervisory Committee MemberWarren Hare, I. K. Barber School Faculty of ScienceSupervisory Committee MemberKlaske Van Heusden, School of EngineeringUniversity ExamineriiAbstractComputing the minima of a convex function is much easier than com-puting the same for a nonconvex function. Every nonconvex function hasthe same minimum as its convex envelope and the minima of the nonconvexfunction are contained in the set of minima of its convex envelope. We com-pute the conjugate of a bivariate piecewise linear-quadratic (PLQ) functionas a first step toward the computation of the convex envelope.Our algorithm starts with computing the convex envelope of each pieceobtaining a rational function defined over a polyhedral subdivision. Thenwe compute the conjugate of each of those pieces and obtain a fractionalform defined over a parabolic subdivision. The last step is to compute themaximum of all those functions to obtain the conjugate of the original piece-wise linear-quadratic function as a piecewise function defined on a parabolicsubdivision.Our implementation in MATLAB uses symbolic computation and ratio-nal numbers to avoid any floating-point errors.We have developed the theory to compute the conjugate by computingthe maximum of the conjugates corresponding to all the pieces and haveimplemented the code to compute the convex envelope, conjugate of eachpiece and the conjugate of the PLQ function.iiiLay SummaryThis research gives an algorithm to compute the conjugate of a bivariatepiecewise linear-quadratic (PLQ) function. This is done is three steps, westart with computing the convex envelope of each piece obtaining a rationalfunction defined over a polyhedral subdivision. Then we compute the con-jugate of each of those pieces and obtain a fractional form defined over aparabolic subdivision. The last step is to compute the maximum of all thosefunctions to obtain the conjugate of the original piecewise linear-quadraticfunction as a piecewise function defined on a parabolic subdivision. We haveimplemented a code in MATLAB using symbolic computation and rationalnumbers to compute the conjugate. The last step to compute the maximumof the conjugates of each division to obtain the conjugate of the originalpiecewise linear-quadratic function is the main focus of this research.ivPrefaceThis thesis is original, unpublished work by the author, Tanmaya Kar-markar conducted under the guidance of Dr. Yves Lucet. No aspects of thiswork involved the use of Generative AI. All images included in this thesisare created by the author.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xivDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3Chapter 2: Preliminaries . . . . . . . . . . . . . . . . . . . . . . 5Chapter 3: Algorithm to compute the conjugate of a PLQfunction. . . . . . . . . . . . . . . . . . . . . . . . . . 123.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.2 Algorithmic approach . . . . . . . . . . . . . . . . . . . . . . 16Chapter 4: Code architecture . . . . . . . . . . . . . . . . . . . 174.1 Data structures and classes . . . . . . . . . . . . . . . . . . . 174.1.1 Classes for piecewise functions . . . . . . . . . . . . . 17viTABLE OF CONTENTS4.1.2 Classes specific to obtaining the conjugate of a PLQfunction . . . . . . . . . . . . . . . . . . . . . . . . . . 224.2 Some methods from class region . . . . . . . . . . . . . . . . . 234.2.1 Intersection of two regions . . . . . . . . . . . . . . . . 234.2.2 Simplifying the region . . . . . . . . . . . . . . . . . . 234.2.3 Set difference (minus) . . . . . . . . . . . . . . . . . . 334.2.4 Maximum of two convex functions over a region . . . 33Chapter 5: Convex envelope of each piece . . . . . . . . . . . . 385.1 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 425.1.1 Subproblems . . . . . . . . . . . . . . . . . . . . . . . 42Chapter 6: Conjugate of each piece . . . . . . . . . . . . . . . . 476.1 Domain of the conjugate . . . . . . . . . . . . . . . . . . . . . 476.1.1 Interior . . . . . . . . . . . . . . . . . . . . . . . . . . 486.1.2 Vertices . . . . . . . . . . . . . . . . . . . . . . . . . . 506.1.3 Edges . . . . . . . . . . . . . . . . . . . . . . . . . . . 536.2 Conjugate Expressions . . . . . . . . . . . . . . . . . . . . . . 556.3 Conjugate - Algorithmic approach . . . . . . . . . . . . . . . 56Chapter 7: Conjugate of a PLQ function . . . . . . . . . . . . . 627.1 Method for Step 3 of the main algorithm. . . . . . . . . . . . 627.2 Step 3a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 647.2.1 Step 3a(i) . . . . . . . . . . . . . . . . . . . . . . . . . 647.2.2 Step 3a(ii) . . . . . . . . . . . . . . . . . . . . . . . . . 697.2.3 Step 3a(iii) . . . . . . . . . . . . . . . . . . . . . . . . 757.2.4 Overall algorithm for Step 3a. . . . . . . . . . . . . . . 757.3 Step 3b . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 767.3.1 Divisions . . . . . . . . . . . . . . . . . . . . . . . . . 79Chapter 8: Experiments . . . . . . . . . . . . . . . . . . . . . . . 838.1 Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 838.1.1 Experiment 1 . . . . . . . . . . . . . . . . . . . . . . . 838.1.2 Experiment 2 . . . . . . . . . . . . . . . . . . . . . . . 838.2 Timing Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . 85Chapter 9: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 909.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 909.1.1 Biconjugate . . . . . . . . . . . . . . . . . . . . . . . . 909.1.2 Step 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . 90viiTABLE OF CONTENTS9.1.3 Step 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . 919.1.4 Parallel Computing . . . . . . . . . . . . . . . . . . . . 919.1.5 Preprocessing . . . . . . . . . . . . . . . . . . . . . . . 919.1.6 Overall Domain . . . . . . . . . . . . . . . . . . . . . . 91Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97Appendix A: Code demonstration . . . . . . . . . . . . . . . . 97Appendix B: Symbolic verification for proofs in Chapter 7. . 99Appendix C: Code - Testing . . . . . . . . . . . . . . . . . . . 102C.1 Unit testing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102C.2 Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102C.3 Overall testing . . . . . . . . . . . . . . . . . . . . . . . . . . 102Appendix D: Derivations for Step 1 . . . . . . . . . . . . . . . 103D.1 Steps to solve the subproblems. . . . . . . . . . . . . . . . . . 103D.1.1 Step (i) . . . . . . . . . . . . . . . . . . . . . . . . . . 103D.1.2 Step (ii) . . . . . . . . . . . . . . . . . . . . . . . . . . 105D.1.3 Step (iii) . . . . . . . . . . . . . . . . . . . . . . . . . 106D.1.4 Step (iv) . . . . . . . . . . . . . . . . . . . . . . . . . . 108viiiList of TablesTable 7.1 Possible divisions of domain obtained in Step 3a (ii). 74Table 8.1 Computational time in seconds: Number of Pieces = 1. 85Table 8.2 Computational time in seconds. Number of Pieces =4; each piece is split in 2 convex pieces. . . . . . . . . 88Table 8.3 Computational time in seconds. Number of Pieces =9; each piece is split in 2 convex pieces. . . . . . . . . 88Table 8.4 Timing summary . . . . . . . . . . . . . . . . . . . . . 88ix.List of FiguresFigure 2.1 Univariate piecewise function from Example 2.4 (inblue) and its convex envelope (in red). . . . . . . . . . 6Figure 2.2 Bivariate PLQ function defined in Example 2.5. . . . 7Figure 3.1 Method to compute the conjugate of a PLQ function. 13Figure 3.2 Illustration of the steps taken on the domains to com-pute the biconjugate of a PLQ function. Note thatthe biconjugate has a polyhedral subdivision, but anunknown explicit formula leading us to label the lastbox “Biconjugate ??”. . . . . . . . . . . . . . . . . . . 14Figure 3.3 Steps to compute the conjugate of a PLQ Functionillustrated on the domain of the respective functions. 15Figure 4.1 Class diagram; arrows represent aggregation relation-ship, i.e., the class where the arrow originates from“contains” the objects of the class the arrow is point-ing to. . . . . . . . . . . . . . . . . . . . . . . . . . . . 18Figure 4.2 Graph of the function r = (155x− 5y+ 4xy+35x2+5y2 − 100)/ (7x− y + 40) . . . . . . . . . . . . . . . . 20Figure 4.3 Region object defined in Example 4.2. . . . . . . . . . 21Figure 4.4 piecewiseFunction object defined in Example 4.3. . . 22Figure 4.5 Intersection of 2 polyhedral sets R = R1 ∩R2. . . . . 25Figure 4.6 Simplify Closed Region (x− 2y = 0 is removed). . . . 26Figure 4.7 Remove inequalities that do not pass through a fea-sible vertex (x+ 7y + 4 = 0 is removed). . . . . . . . 27Figure 4.8 Process of simplifying an unbounded region with morethan two inequalities through a vertex. . . . . . . . . 29Figure 4.9 Simplifying a region with a quadratic inequality: Ex-ample 4.5. . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 4.10 Singleton region. . . . . . . . . . . . . . . . . . . . . . 33Figure 4.11 Minus with one output region, R = R1 −R2 . . . . . 35xLIST OF FIGURESFigure 4.12 Minus with two output regions (R1 −R2 = R3 +R4). 37Figure 5.1 Polyhedral Subdivision. . . . . . . . . . . . . . . . . . 39Figure 5.2 Function restricted to edge. . . . . . . . . . . . . . . . 40Figure 5.3 Bilinear form with domain as defined in Example 5.1. 42Figure 5.4 Polyhedral Subdivision of domain for Example 5.1. . 45Figure 5.5 Convex Envelope for Example 5.1. . . . . . . . . . . . 45Figure 6.1 Step 2: Computing the conjugate of a convex bivari-ate rational function over a polytope. . . . . . . . . . 48Figure 6.2 Conjugate Domain for Example 6.9 . . . . . . . . . . 55Figure 7.1 Step 3: Computing the conjugate of a PLQ function. 63Figure 7.2 Intersection of domains of two conjugates. . . . . . . 65Figure 7.3 Computing dom(si,1,k)∩dom(si,2,k) using dom(ri,1)∩dom(ri,2). . . . . . . . . . . . . . . . . . . . . . . . . 68Figure 7.4 Further division by a parabolic curve as explained inExample 7.4. . . . . . . . . . . . . . . . . . . . . . . . 70Figure 7.5 Merging adjacent pieces with the same function. . . . 75Figure 7.6 Maximum Conjugate for piece 1. . . . . . . . . . . . . 77Figure 7.7 Maximum Conjugate for piece 2. . . . . . . . . . . . . 78Figure 7.8 Maximum Conjugate over all pieces. . . . . . . . . . . 80Figure 7.9 Conjugate of the PLQ function from Example 7.13. . 82Figure 8.1 Experiment 1: Computing the conjugate of a functionwith the domain stored as one piece gives the sameresult as the same function with the domain storedas two pieces. . . . . . . . . . . . . . . . . . . . . . . 84Figure 8.2 Domains used in Experiment 2. . . . . . . . . . . . . 86Figure 8.3 Polyhedral subdivision of domain of convex envelopefor 8.2(c). . . . . . . . . . . . . . . . . . . . . . . . . . 87Figure 8.4 Conjugate of xy over unit square with each piece sub-divided into n2 pieces and corresponding computationtime per piece and per edge. . . . . . . . . . . . . . . 87Figure 9.1 Domain of a conjugate illustrating the potential forimprovement. The conjugates corresponding to ver-tices and edges that are not a part of the divisionwhen computing the convex envelope of each pieceappear exactly the same in the final conjugate. . . . . 92xiLIST OF FIGURESFigure A.1 Function f(x, y) = xy with piecewise polyhedral do-main. . . . . . . . . . . . . . . . . . . . . . . . . . . . 97.xiiList of Algorithms1 Computing the conjugate of a PLQ function . . . . . . . . . . 162 Computing the intersection of two regions: R = R1 ∩R2 . . . 243 Simplifying a closed region . . . . . . . . . . . . . . . . . . . . 264 Removing inequalities that do not pass through a feasible vertex 275 Simplifying an unbounded region . . . . . . . . . . . . . . . . . 326 Computing the set difference between two regions . . . . . . . 347 Computing the maximum of two functions over a region . . . . 368 Computing the parabolic curve from rational expressions . . . 499 Computing the normal Cone at each vertex . . . . . . . . . . . 5010 Computing the subdifferential at vertices where limits exist . . 5111 Computing the subdifferential at vertices where limits do notexist (by continuity) . . . . . . . . . . . . . . . . . . . . . . . . 5212 Computing the normal cone at points on edges. . . . . . . . . 5313 Computing the subdifferential at points on edges with verticesof 00 type . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5414 Computing the conjugate expression at vertices . . . . . . . . . 5615 Computing the quadratic expression at edges when the coeffi-cient of y is 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . 5716 Computing quadratic expressions at edges where ξ21 +mξ21 = 0 5817 Computing the fractional expression at edges with ξ21+mξ21 ̸= 0 5918 Computing the conjugate expression at edges . . . . . . . . . . 6019 Computing the conjugate of rational expressions . . . . . . . . 6120 Computing the Cartesian product of two lists of regions . . . . 6621 Computing the intersection of two lists of dual regions that areadjacent in the primal domain . . . . . . . . . . . . . . . . . . 6722 Computing the maximum of two functions over a given domain 7123 Computing the maximum to obtain the conjugate of a piece . 7624 Computing the maximum to obtain the conjugate of the PLQfunction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79xiiiAcknowledgementsI extend my deepest gratitude to my supervisor, Dr. Yves Lucet, forentrusting me with this project. His meticulous attention to detail and gen-erous mentorship have been invaluable throughout my research and writingprocess. His guidance has not only shaped my work but also my growth asa researcher.I am immensely thankful to the members of my supervising committee,Dr. Heinz Bauschke and Dr. Warren Hare, for their insightful feedbackon my thesis. Under Dr. Bauschke’s guidance, exploring convex analysishas been both enriching and integral to this project. Dr. Hare’s expertisein Derivative Free Optimization provided a fresh perspective on algorithmsand convergence analysis, which greatly enhanced the depth of my research.I am also thankful to all of my committee members for making my de-fense a memorable milestone. I would also like to mention the externalexaminer, Klaske Van Heusden for her insightful questions that would defi-nitely help in continuing this research.My gratitude also extends to my husband, Dr. Narendra Karmarkar,whose mentorship has been instrumental in honing my skills in mathematicalsoftware development. Additionally, I am indebted to my seniors, Shamb-havi Singh, Zenab Kagdiwala, Namrata Kundu, and Sayan Sadhukhan, fortheir unwavering support and guidance. A special mention goes to my col-league, Meghana Khuntia, for her constant support over the past two years.Lastly, I am grateful to my lab mates, Saira, Paula, Paavai, and Pragya,whose camaraderie made the journey all the more enjoyable.Lastly, I would like to thank The University of British Columbia, forfinancial assistance. This research was supported by Natural Sciences andEngineering Research Council of Canada (NSERC), the University of BritishColumbia, Okanagan and was performed in the Computer-Aided ConvexAnalysis (CCA) library funded by a Leaders Opportunity Fund (LOF, JohnR. Evans Leaders Fund — Funding for research infrastructure) from theCanada Foundation for Innovation (CFI) and by British Columbia Knowl-edge Development Fund (BCKDF). Their consistent financial support al-lowed me to dedicate myself fully to this research endeavor.xivTo Ame´y and Siddhant.xvChapter 1IntroductionMathematical optimization is the selection of a best element, with regardto some criterion, from some set of available alternatives. More generallya solution to an optimization problem consists of a process or algorithmto maximize or minimize a real valued function by systematically choosingpoints from within an allowed set and computing the value of the function.The generalization of optimization theory and techniques constitute a largearea of applied mathematics. Optimization problems arise in all quantitativedisciplines from computer science and engineering to operations researchand economics. The development of solution methods has been of interestin mathematics for centuries.1.1 Literature reviewGlobal optimization is a branch of applied mathematics and numericalanalysis that attempts to find the global minima or maxima of a functionor a set of functions on a given set. Global optimization is distinguishedfrom local optimization by its focus on finding the minimum or maximumover the given set, as opposed to finding local minima or maxima. Findingan arbitrary local minimum is relatively straightforward by using classicallocal optimization methods. Finding the global minimum of a function isfar more difficult; analytical methods are frequently not applicable, and theuse of numerical solution strategies often leads to very hard challenges.A large number of algorithms proposed for solving the nonconvex prob-lems are not capable of making a distinction between locally optimal solu-tions and globally optimal solutions, and treat the former as actual solutionsto the original problem. Global optimization is the branch of applied math-ematics and numerical analysis that is concerned with the development ofalgorithms that are capable of guaranteeing convergence in finite time to theactual optimal solution of a nonconvex problem.Convex optimization is a subfield of mathematical optimization thatstudies the problem of minimizing convex functions over convex sets [NN94].Many classes of convex optimization problems admit polynomial-time algo-11.1. LITERATURE REVIEWrithms, whereas mathematical optimization is in general NP-hard [MK87;Sah74].In many practical problems the function to be minimized is noncon-vex and this is difficult. Computing the (global) minima of a convex func-tion [BC17, Theorem 16.2] is much easier than computing the same for anonconvex function. However every nonconvex function shares the sameminimum as its convex envelope.Computing the convex envelope of a function is a hard problem; evencomputing the convex envelope of a multilinear function over a unit hyper-cube is NP-Hard [Cra89]. However, results for specific functions exist in theliterature, in particular for quadratic bivariate polynomials. The problem ofcomputing the convex envelope for any quadratic bivariate polynomial canbe transformed into computing the convex envelope of the sum of a bilinearpolynomial and a linear function. Convex envelopes for bilinear functionsover triangles, rectangles and special polytopes have been studied in [AF83;SA90; Lin05; AB10; Ans12].We work with a class of functions called Piecewise Linear-Quadraticfunctions. These are piecewise quadratic functions defined on a union ofpolyhedral sets. They are abbreviated as PLQ functions.Since the closed convex envelope of a function is equal to its bicon-jugate, computing the conjugate provides another method to obtain theconvex envelope. The early idea of the computation of convex transformscan be traced back to [Mor65]. However, development of most of the al-gorithms in computational convex analysis began with the Fast LegendreTransform (FLT) [Bre89], which is studied in [Cor96; Luc96]. A linear-time algorithm is introduced in [Luc97]. Extension to nongrid domains areproposed in [HL18; GL13] with a generalization to the partial conjugatein [GKL14]. Computation of the conjugate of convex univariate PLQ func-tions has been well studied and linear time algorithms have been developed,both for the PLQ [LBT09] and the GPH models [GL11].Implementation of conjugate computation algorithms can be found in theCCA library [Luc17; Luc21] that contains efficient algorithms to manipulateconvex functions and to compute fundamental convex analysis transformsarising in the field of convex analysis. Recent algorithms to compute theconjugate numerically have been based on either parameterization [HL07],manipulation of graphs (GPH model) [GL13], or the computation of theMoreau envelope [Luc06]. More complex operators such as the proximalaverage operator [Bau+08; BLT08] can be built by using a combination ofaddition, scalar multiplication,and conjugacy operations.Computational Convex Analysis has many applications in the fields of21.2. MOTIVATIONimage processing, network communication, PDE, geographic informationsystems, computer aided design, molecular biology, medical imaging, com-puter graphics and robotics; see the survey [Luc10].PLQ functions play a significant role in variational analysis [RW98, Sec-tion 10E, Section 11D, p. 440] due to the availability of calculus rules [RW98,Example 11.28, Proposition 11.32, Corollary 11.33, Proposition 12.30], andtheir duality property [RW98, Theorem 11.42, Example 11.43, Theorem11.42]. Compared to piecewise linear functions, PLQ functions capture thecurve of the original function more accurately and are continuous but maynot be continuously differentiable. On each piece they are second-order dif-ferentiable in the interior. The set of PLQ functions is closed under commonconvex operators; see [RW98, Page 484] and [GKL14].1.2 MotivationWhile the convex envelope of a piecewise function is not always the sameas the convex envelope of each piece, for bivariate quadratic polynomials wecan leverage the convex envelope of each piece to obtain the convex envelope.Let us denote the ith piece of function f by (fi, Pi) where fi is a function andPi is a set. The indicator function of C ⊂ Rn is denoted IC : R→ (−∞,∞]with IC(x) = 0 if x ∈ C and +∞ otherwise. The Fenchel conjugate of afunction f : Rn → R ∪ {∞} is denotedf∗(s) = supx∈Rn[sTx− f(x)].(It is also called the convex conjugate, or simply conjugate.) The Legendre-Fenchel transform is the operator that transforms f 7→ f∗, whereas theconjugate f∗ : s 7→ f∗(s) is the function obtained by that transform. For afunction f , the pointwise supremum of all convex underestimators of f , iscalled the Convex Envelope and denoted by conv(f).From [HL93, Theorem 3.4.1], we have (infi fi)∗ = supi f∗i . Using [HL93,Proposition 2.6.1] and the biconjugate theorem conv f = f∗∗, we getconv[infi(fi + IPi)] = conv[infi(conv(fi + IPi))],=[infi(conv(fi + IPi))]∗∗,=[supi([conv(fi + IPi)]∗)]∗.Hence we can find the conjugate of a bivariate PLQ function by followingthe following three steps:31.2. MOTIVATION1. Compute the convex envelope of each piece conv(fi + IPi).2. Compute the conjugate of each convex piece [conv(fi + IPi)]∗.3. Compute the maximum of the conjugates over the entire PLQ functionto obtain f∗ = supi(conv(fi + IPi)∗).Step 1 of this process, finding the convex envelope of each piece, is com-puted as explained in [Loc16] where a method to compute the convex en-velopes of bilinear forms over general polytopes is given (using convex anal-ysis, computing the convex envelope of a general quadratic function reducesto computing the convex envelope of the bilinear function (x, y) 7→ xy).Kumar [Kum19] showed how to compute the conjugate of a rationalfunction defined on a polytope, which is Step 2 of the above process.In this thesis, we explain how to perform Step 3 by computing the maxi-mum of the functions obtained in Step 2 thereby obtaining the conjugate of abivariate PLQ function. Repeating steps 2 and 3 to compute the biconjugatef∗∗ = (supi(conv(fi + IPi)∗))∗ is left for future work.Our implementation in MATLAB uses symbolic computation and ratio-nal numbers to avoid any floating-point errors. Using floating-point arith-metic can give rise to partitions of the domain with degenerate subsets, andfloating-point coefficients of functions make it difficult to detect equal func-tions on adjacent domains. After observing these difficulties, we decided towork with symbolic representation of our functions with rational coefficients.We begin by giving some preliminaries in Chapter 2 and the overallalgorithm in Chapter 3. Chapter 4 discusses the implementation followedby details of the first three steps of the algorithm, in chapters 5, 6 and 7.Chapter 8 contains some experimental results followed by the conclusionand future work in Chapter 9.4Chapter 2PreliminariesIn this chapter, we recall known results, define terms used in the rest ofthe thesis, and setup our notations.First, we want to avoid manipulating degenerate functions like f(x) =∞for all x, or a function for which there is x ∈ Rn with f(x) = −∞.Definition 2.1. A function f : Rn → R ∪ {+∞} is called proper if thereexists x ∈ Rn with f(x) <∞.Now we define the class of functions we work with.Definition 2.2 (Piecewise functions). A function f : Rn → R ∪ {∞} iscalled a piecewise function if it can be written asf(x) =f1(x), x ∈ R1,f2(x), x ∈ R2,...fn(x), x ∈ Rn.where Ri is the domain of function fi with Ri ∩ Rj = ∅ for i ̸= j, and∪iRi = Rn.Definition 2.3 (Piece). We call each component of the piecewise function,(fi, Ri) a piece of the piecewise function.A simple example of a nonconvex piecewise function is provided next.Example 2.4. The univariate nonconvex piecewise function,f(x) =f1(x) = x2 if x ≤ 0,f2(x) = 1− (x− 1)2 if 0 < x < 1,f3(x) = x2 if x ≥ 1.is illustrated in Figure 2.1. The blue curve represents the second piece, whilethe first and third piece are represented by the red curve.5Figure 2.1: Univariate piecewise function from Example 2.4 (in blue) andits convex envelope (in red).In this thesis, we focus on bivariate functions similar to the next example.Example 2.5. The bivariate nonconvex piecewise function,f(x) ={f1(x, y) = x2 + xy − y2 if (x, y) ∈ P1,f2(x, y) = 2x2 − xy − y2 if (x, y) ∈ P2,whereP1 = {(x, y) : −x+ 2y ≤ 0 and 3x− 2y − 12 ≤ 0 and −x− 2y + 4 ≤ 0},P2 = {(x, y) : x− 2y ≤ 0 and − 4x+ y + 7 ≤ 0 and 2x+ 3y − 21 ≤ 0},is shown in Figure 2.2.Definition 2.6 ([RW98, Definition 10.20] Piecewise Linear-Quadratic Func-tions). In Definition 2.2, if all functions, f1, f2, . . . , fn are quadratic (or lin-ear) and all the sets Ri are polyhedral, then the function is called a PiecewiseLinear-Quadratic function (PLQ function).Example 2.5 is an example of a PLQ function.We now recall the definition of convex sets and convex functions.Definition 2.7 (Convex Set). A set S ⊂ Rn is convex if∀x, y ∈ S,∀λ ∈ [0, 1], λx+ (1− λ)y ∈ S.6Figure 2.2: Bivariate PLQ function defined in Example 2.5.Definition 2.8 (Convex Function). A function f : Rn → R∪{∞} is convexif and only if dom f = {x : f(x) ∈ R} is a convex set, and for all x, y ∈ dom f ,λ ∈ [0, 1], we havef(λx+ (1− λ)y) ≤ λf(x) + (1− λ)f(y).It is strictly convex if the inequality is strict for all λ ∈ (0, 1).Next we define terms leading to the convex envelope and its conjugate.Definition 2.9 (Lower Semicontinuous Functions). Let f : Rd → [−∞,∞].Then f is lower semicontinuous (lsc) at x ∈ Rd if for every sequence (xn)n∈N ∈Rd, we havelimn∈Nxn = x =⇒ f(x) ≤ lim infnf(xn).Definition 2.10 (Convex underestimator). Any convex function that liesbelow a given function is called an underestimator.Definition 2.11 (Closed Convex envelope). The closed convex envelope isdefined as the largest lsc convex underestimator of f .Example 2.12. The convex envelope of the univariate nonconvex piecewisefunction in Example 2.4 isconv(f(x)) =f1(x) = x2 if x ≤ 0,f2(x) = x if 0 < x < 1,f3(x) = x2 if x ≥ 1.(2.1)It is illustrated by the red curves in Figure 2.1.7Definition 2.13 (Fenchel Conjugate). Let f : Rn → [−∞,+∞] Thenf∗ : Rn → [−∞,+∞] : y 7→ supx∈Rn(xT y − f(x))is the Fenchel conjugate of f .Definition 2.14 (Biconjugate). The biconjugate is defined as (f∗)∗ andnoted as f∗∗.Fact 2.15. [RW98, Theorem 11.1] When f is proper, closed and convex,f = f∗∗, otherwise f∗∗ is the closed convex envelope of f .We now define terms that are required to describe the domains of thepieces of our PLQ functions.Definition 2.16 (Hyperplanes). A hyperplane H ⊂ Rd is a set of the formH = {x ∈ Rd : αTx − β = 0}, where the vector α ∈ Rd\{0} is called anormal vector to H.Definition 2.17 (Halfspaces). A halfspace H ⊂ Rd is a set of the formH = {x ∈ Rd : αTx − β ≤ 0} or H = {x ∈ Rd : αTx − β ≥ 0} where thevector α ∈ Rd\{0} is called a normal vector to H.Definition 2.18 (Polytope). A polytope is a convex combination of finitenumbers of vertices V = {v1, . . . , vk} and it is expressed as,conv V ={λ1v1 + · · ·+ λkvk : λi ≥ 0,k∑i=1λi = 1}={(v1, . . . , vk)λ : λi ≥ 0,k∑i=1λi = 1}={V λ : λi ≥ 0,k∑i=1λi = 1}.where V is the matrix formed by column vectors.Note that by our definition, polytopes are always bounded.Definition 2.19 (Polyhedral cone). A polyhedral cone is defined as a coniccombination of a finite numbers of directions C = {d1, . . . , dm}; we denotecone C = {µ1d1 + · · ·+ µmdm : µi ≥ 0}= {(d1, . . . , dm)µ : µi ≥ 0}= {Dµ : µi ≥ 0}.where D is the matrix formed by column vectors.8Definition 2.20 (Polyhedral set). A polyhedral set or polyhedron P ⊂ Rdis described as the intersection of a finite number of closed halfspaces andhyperplanes; it has the H-representation P = {x ∈ Rd : Ax ≤ b} withA ∈ Rm×d and b ∈ Rm.Fact 2.21. We can also describe a polyhedron in Vertex-representation. Ac-cording to the Minkowski-Weyl theorem [RW98, Corollary 3.53], a polyhe-dron can be written as a sum of a polytope and a polyhedral cone. Therefore,in V -representation a polyhedron P is described asP = conv V + coneC={V λ+Dµ : λ ∈ Rk, λ ≥ 0, µ ≥ 0,k∑i=1λi = 1}.Fact 2.22. Since the intersection of convex sets is convex, polyhedral setsare convex.We recall that the dimension of a set S ⊂ Rd is the dimension of thesmallest affine subspace containing S.Definition 2.23 (Face of a polyhedron [Roc70, Page 162]). A d-dimensionalface of a convex set C ⊂ Rn is a d-dimensional convex subset C ′ of C suchthat every (closed) line segment in C with a relative interior point in C ′has both endpoints in C ′. The empty set and C itself are faces of C. Thetwo-dimensional faces of C are called faces.Definition 2.24 (Edge of a polyhedron [Roc70, Page 162]). The one-dimensional faces of a convex set C are called edges. For u, v ∈ Rd, u ̸= vand δ = (δ1, δ2) ∈ ∆ ⊂ R2, an edge can be written as{x ∈ Rd : x = δ1u+ δ2v, δ1 + δ2 = 1}.An edge in R2 is a1. Segment if ∆ = R+ R+ where R+ = {α ∈ R : α ≥ 0},2. Ray if ∆ = R+ R,3. Line if ∆ = R2.Definition 2.25 (Vertex of a polyhedron ). [Roc70, Page 162] The zero-dimensional faces of a convex set C are called vertices and are the extremepoints of C.9Definition 2.26 (Polyhedral Subdivision). A convex set defined as theunion of a finite number of polyhedral regions, namelyR =⋃i∈{1,2,..n}Ri,where R ⊆ R2, is said to be a polyhedral subdivision, if Ri is a polyhedralset and for any j, k ∈ {1, . . . , n}, j ̸= k,Rj ∩Rk is either empty, a vertex oran edge.Polyhedral subdivisions eliminate degenerate subdivisions for which theintersection of two polyhedral regions is a strict subset of an edge, or of aface.Definition 2.27 (Entity). Assume that f : R2 → R ∪ {+∞} is a PLQfunction and that dom f = ∪iPi is a polyhedral subdivision. An entity is ad-dimensional face of Pi for some index i. An entity is either a vertex, anedge or a face of the polyhedral subdivision of dom f [GL13].We now define the geometric shapes required to describe the domain ofthe conjugate functions we obtain.Definition 2.28 (Parabola, Parabolic arc/curve). A parabola P ⊂ R2 is asubset of the plane that can be written asP = {(x, y) ∈ R2 : ax2 + bxy + cy2 + dx+ ey + f = 0},where a, b, c, d, e, f ∈ R are not all zero and satisfy b2 − 4ac = 0 .Note that our definition includes lines when a = b = c = 0, and theempty set when a = b = c = d = e = 0, f ̸= 0; but excludes the entire planesince we impose that (a, b, c, d, e, f) ̸= 0.Definition 2.29 (Parabolic Region). A parabolic region Pr ⊂ R2 is formedby the intersection of a finite number of parabolic inequalities. It can bewritten asPr = {x ∈ R2 : aix2 + bixy + ciy2 + dix+ eiy + fi ≤ 0, i = 1, . . . , k},where (ai, bi, ci, di, ei, fi) ̸= 0, and b2i − 4aici = 0.The following sets are special cases of parabolic regions: polyhedral sets,polyhedral cones, polytopes, edges, and vertices. The set {(x, y) ∈ R2 : y ≥x2} is an example of a non-polyhedral parabolic region.10Definition 2.30 (Parabolic subdivision). A convex set R ⊂ R2 is called aparabolic subdivision if R can be written as the finite union of parabolicregions Ri, i.e., R = ∪i∈{1,2,..n}Ri and for any j, k ∈ {1, . . . , n}, j ̸= k, theintersection Rj ∩Rk is either empty or is contained in a parabola.A polyhedral subdivision is a special case of a parabolic subdivision(since lines are a special case of parabolas).Definition 2.31 (Face of a parabolic region). A face of a parabolic regionPr ⊂ R2 is defined as the interior of a nonempty intersection of a finitenumber of parabolic inequalities. It can be written asPr = {x ∈ R2 : aix2 + bixy + ciy2 + dix+ eiy + fi < 0, i = 1, . . . , k},where (ai, bi, ci, di, ei, fi) ̸= 0, and b2i − 4aici = 0.Finally, we recall the definition of subdifferentials and normal cones sincethey are required in the computation of the conjugate. We refer to [BM23]for the notations and similar concepts.Definition 2.32 (Subdifferential). The subdifferential of a function f atany x ∈ dom(f) is∂f(x) = {s : f(y) ≥ f(x) + sT (y − x), ∀y ∈ dom f}.When x /∈ dom(f), we define ∂f(x) = ∅.We call subgradient any point s in ∂f(x).Definition 2.33 ([Roc70, p.15] Normal cone). The normal cone at x ∈ Rnof a nonempty subset C ⊂ Rn isNC(x) ={{u ∈ Rn : (∀c ∈ C), (c− x)Tu ≤ 0} if x ∈ C;∅ if x ̸∈ C. (2.2)11Chapter 3Algorithm to compute theconjugate of a PLQ function.In this chapter we give the algorithm to compute the conjugate of thePLQ function. We also define the notation and give some illustrative exam-ples explaining the approach we are taking in computing the conjugate.3.1 NotationWe use the symbols q for quadratic functions, r for rational functionsand s for fractional forms. We now fix our notations for each intermediatecomputation by naming each of the piecewise function involved.The original PLQ function and its domain are noted f and dom f wheref = mini=1...nq(qi + Idom qi).f has nq pieces where qi is the quadratic function defining the ith piece.The convex envelope of qi + Idom qi is denotedconv(qi + Idom qi) = minj=1...nri(ri,j + Idom ri,j ).Here ri,j represents a rational function that defines the jth piece of theconvex envelope of the ith piece of f .The conjugate of ri,j + Idom ri,j is written as(ri,j + Idom ri,j )∗ = mink=1,...nsi,j(si,j,k + Idom si,j,k).si,j,k represents a fractional function which is the kth piece of r∗i,j .The conjugate of conv(qi + Idom qi) is[conv(qi + Idom qi)]∗ = maxj=1...nrimink=1,...nsi,j(si,j,k + Idom si,j,k),= mink=1...nsi(si,k + Idom si,k).123.1. NOTATIONFigure 3.1: Method to compute the conjugate of a PLQ function.The conjugate of f is denotedf∗ = maxi=1...nq(si,k + Idom si,k) = mink=1...ns(sk + Idom sk).We illustrate the overall flow of computing the conjugate of the PLQfunction in Figures 3.1 and 3.2. We are given a PLQ function with n pieces.Each piece has a nonconvex bivariate function defined over a polyhedralregion. As an example in Figure 3.3(a), we have a PLQ function with twopieces.The algorithm follows the following stepsStep 1: Given qi compute the convex envelope of qi+ Idom qi to obtain ri,j .Step 2: From ri,j , compute the conjugate of ri,j + Idom ri,j to obtain si,j,k.Step 3: From si,j,k, compute the conjugate of conv(qi + Idom qi) to obtainsi,k, and deduce sk to obtain the conjugate of f .We apply Step 1 and get the convex envelope for each piece, as shownin Figure 3.3(b) and Figure 3.3(c). This is followed by Step 2 to obtain theconjugate of each convex piece. The conjugates are shown in Figures 3.3(d-g) and finally the maximum of the conjugates are computed in Steps 3.Figures 3.3(h-i) show the conjugate of each piece and finally Figure 3.3(j)shows the conjugate over the entire PLQ function.133.1. NOTATIONFigure 3.2: Illustration of the steps taken on the domains to compute thebiconjugate of a PLQ function. Note that the biconjugate has a polyhedralsubdivision, but an unknown explicit formula leading us to label the lastbox “Biconjugate ??”.143.1. NOTATION(a) qi. (b) r1,j . (c) r2,j .(d) s1,1,k. (e) s1,2,k. (f) s2,1,k. (g) s2,2,k.(h) s1,k. (i) s2,k.(j) sk.Figure 3.3: Steps to compute the conjugate of a PLQ Function illustratedon the domain of the respective functions.153.2. ALGORITHMIC APPROACH3.2 Algorithmic approachTo compute the conjugate of a PLQ function, we follow the first threesteps stated earlier. Algorithm 1 is used to find the conjugate. We iterateover each piece and compute its convex envelope, conjugate and maximum(conjugate of each piece). Finally we compute the maximum over the con-jugates of all the pieces to obtain the conjugate of the entire PLQ function.A Step 4 that calculates the biconjugate has not yet been implementedbut will be added as the last step in the future.Algorithm 1: Computing the conjugate of a PLQ functionInput: p (PLQ function)Output: p∗foreach pi in p do% pi is the ith piece of pqi = pi.convexEnvelope %qi is a piecewise function.foreach ri,j in qi dosi,j = ri,j .conjugate %si,j is a piecewise function.endsi = si,j .maximumConjugate (maximum over j).endp∗ = si.maximumConjugate (maximum over i).16Chapter 4Code architectureConsidering the overall flow as illustrated in Figure 3.3 and the commonfunctionality between the different steps, we have defined numerous classesto organize the code. The overall idea behind structuring this code is to keepthe code as close to the mathematical formulae and algorithms as possible.This would make it easier to use, modify and add. The code is implementedin MATLAB. All computations are done symbolically.Symbolic computing keeps the code close to the mathematical formulaeand hence it is easier to interpret, implement and debug. Although symboliccomputing is significantly slower, as we are working with only two variablesupto degree two, we can use this technique.In MATLAB, Symbolic Math Toolbox introduces a special data type —symbolic objects. This data type includes symbolic numbers, symbolic vari-ables, symbolic expressions, and symbolic functions. It also includes vectors,matrices, and multidimensional arrays of symbolic numbers, variables, ex-pressions, and functions. Using symbolic objects in computations indicatesthat MATLAB must perform these computations analytically instead of nu-merically. Symbolic computations are exact and are not prone to round-offerrors.The code we have implemented is available at“git@github.com:tanmaya11/convex.git”.4.1 Data structures and classesThe main classes used are described in this section. We have othersupporting classes defined in the actual code.4.1.1 Classes for piecewise functionsAs we are dealing with piecewise functions, we first define classes to storeand operate on these functions. The basic class diagram depicting the aboveclasses is given in Figure 4.1.174.1. DATA STRUCTURES AND CLASSESFigure 4.1: Class diagram; arrows represent aggregation relationship, i.e.,the class where the arrow originates from “contains” the objects of the classthe arrow is pointing to.184.1. DATA STRUCTURES AND CLASSESsymbolicFunction classAt the lowest level, we have defined a class called symbolicFunction. Thisclass stores functions symbolically; using a MATLAB symbolic variable.Along with the function, we also store the number of variables and thesymbolic variables used to define the function. We define a class to storea function instead of using MATLAB symbolic variable to store a functiondirectly, in order to handle special cases, (e.g., rational functions evaluatedto 00). This also allows us to isolate methods required at the lowest leveland add more functionality as required. These can be further optimized asrequired without disturbing the rest of the code.Example 4.1. Consider the functionr(x, y) =155x− 5y + 4xy + 35x2 + 5y2 − 1007x− y + 40 . (4.1)Its graph is plotted in Figure 4.2. We first declare symbolic variables x andy in MATLAB and then define r in terms of x and y as given below.x = sym( ’ x ’ ) ;y = sym( ’ y ’ ) ;r = symbol icFunct ion((155x− 5y + 4xy + 35x2 + 5y2 − 100)/(7x− y + 40)) ;This creates a variable r as an object of class symbolicFunction. Inter-nally variable r storesr.nv = 2r.vars = [x, y]r.f = ((155x− 5y + 4xy + 35x2 + 5y2 − 100)/(7x− y + 40))Of these variables nv and vars are private, f currently has public accessbut will be made private in the future. Extracting and storing nv and varsis useful in methods like gradient and Hessian where g = r.gradient wouldreturn an array of size nv with derivatives of r with respect to each variablein vars.region classThis class is used to store a feasible, nonempty set defined by polynomialinequalities. These polynomials are stored as objects of class symbolicFunc-tion. We store a list of inequalities as LHS ≤ 0 where only LHS is stored in194.1. DATA STRUCTURES AND CLASSESFigure 4.2: Graph of the function r = (155x−5y+ 4xy+35x2+5y2−100)/(7x− y + 40)the list. The inequalities are normalized prior to storing to enable compar-isons. We also store the vertices defining the feasible region. We can createan object of type region either by supplying inequalities defining the region,or by supplying vertices in case of a bounded region. When we create anobject of type region, if vertices are given we compute the inequalities and ifinequalities are given we compute the vertices. We are assuming that thereis at least one vertex in the region. If the number of vertices is zero we callit an empty region and return an empty object.Example 4.2. Number of Vertices = 2.Variables: s1, s2.Vertices: v = (13944 ,−4544), (5411 , 811).Intersection of the following inequalities:−s1 − 7s2 − 4 ≤ 0,s1 + 7s2 − 10 ≤ 0,s21 + 14s1s2 + 148s1 + 49s22 − 196s2 − 684 ≤ 0.Figure 4.3 shows the feasible region in blue along with defining inequal-ities.piecewiseFunction classThis class is the representation of a piecewise function. In a piecewisefunction, each piece has a function defined on a specific domain. Hence we204.1. DATA STRUCTURES AND CLASSES.Figure 4.3: Region object defined in Example 4.2.define a class to encapsulate a function and its domain together. We createtwo parallel lists, one for the functions (symbolicFunction) and the other forthe domains (region). We encapsulate these lists together in this class tostore all the pieces of a piecewise function.Example 4.3. This example shows the first piece of a piecewise function.Function: −4x− 5y − 20.Domain:Number of Vertices = 3.Variables: x, y.Vertices: v = (-5,-4),(0,-4),(-5,5).Intersection of the following inequalities:−y − 4 ≤ 0, (4.2)−x− 5 ≤ 0, (4.3)x+59y +209≤ 0. (4.4)This example is illustrated in Figure 4.4. In the illustration, the functionshifted for clarity is shown in green, while the domain is shown in blue.214.1. DATA STRUCTURES AND CLASSESFigure 4.4: piecewiseFunction object defined in Example 4.3.4.1.2 Classes specific to obtaining the conjugate of a PLQfunctionIn this section, we describe the classes we have used to store and operateon a PLQ function in order to compute its conjugate.plq 1piece classFrom this class onwards, we are being more specific to the problem weare solving. This class stores one piece of the PLQ function. It has thefollowing fields:1. Original function and domain.2. Convex envelope - piecewiseFunction.3. Conjugate - piecewiseFunction.4. iaConjugate - indexing to show break up of conjugate list to correspondto each division.5. Max conjugate - piecewiseFunction.This class contains methods defined to find the convex envelope, conju-gate, maximum (conjugate) and other supporting functions for one piece.224.2. SOME METHODS FROM CLASS REGIONplq classThis class stores a PLQ function. It has a list of pieces (objects ofplq 1piece) and the conjugate (piecewiseFunction). The biconjugate (piece-wiseFunction) will also be stored here.4.2 Some methods from class regionSome core operations used in the next three chapters to manipulate thedomain are given here.4.2.1 Intersection of two regionsWe find the intersection of two objects of class region. The methodreturns a feasible intersection or an empty region, when the output is notfeasible; Algorithm 2 lists the steps while an example is shown in Figure 4.5.In this algorithm we return the intersection when the interior of the inter-section is nonempty. In order to compute the intersection R1 ∩R2, we firstput all the inequalities defining R1 and R2 together in a list. If there areany inequalities; f ≤ 0 and −f ≤ 0, the intersection would give f = 0.As we do not store degenerate regions, if we encounter degenerate regionswe return an empty region. If there are no degenerate regions; we invokethe constructor for region with the list of inequalities. The constructor willreturn the region empty if there are no feasible vertices else it returns anobject of class region.4.2.2 Simplifying the regionAs seen in the output of Figure 4.5, when we intersect two regions, we canhave a region which has inequalities that are not required (like x+ 4 ≤ 0 inFigure 4.5(c)) and hence can be removed. We call removing these redundantinequalities simplifying the region. We have currently defined two differentalgorithms for simplifying regions. The first algorithm works only when theregion is bounded. The second one works with unbounded regions.Simplify closed region.In Step 1, all the regions are closed. When we create any region, wealways compute all its vertices. In order to remove redundant inequalities,we remove one linear inequality at a time and check that all the vertices234.2. SOME METHODS FROM CLASS REGIONAlgorithm 2: Computing the intersection of two regions: R =R1 ∩R2Input: R1, R2Output: R% Copy inequalities from R1 into a list ln = 0foreach ineq ∈ R1 don = n+ 1l(n) = ineqendforeach ineq ∈ R2 do% Check if an inequality in R2 contradicts one in R1if −ineq ∈ R1 thenR = ∅returnend% Do not store duplicate inequalitiesif ineq ∈ R1 thencontinueend% Add inequalities of R2 to ln = n+ 1l(n) = ineqend% Create a region with inequalities in list l.R = region(l)% If region is infeasible - constructor of region returns empty.244.2. SOME METHODS FROM CLASS REGION 7 4 095−5 ≤(a) R1.y+5 0+40(b) R2.(c) R.Figure 4.5: Intersection of 2 polyhedral sets R = R1 ∩R2.remain the same. Algorithm 3 shows the pseudo-code, which is illustratedin Figure 4.6. This algorithm is linear time in the number of inequalities.Simplify unbounded regionIn Steps 2 and 3, our domains are unbounded regions. Here we also havequadratic inequalities. In order to simplify these regions, we work in twosteps. We also assume that we have at least one vertex, as that is alwaysthe case in the problem we are solving (in addition, one can always add anarbitrary vertex to a region if needed by subdividing the region).Simplify region with inequalities not passing through feasiblevertices.We want to remove the inequalities whose removal will not change thefeasible region. In Algorithm 4 we remove all the inequalities that do not254.2. SOME METHODS FROM CLASS REGIONAlgorithm 3: Simplifying a closed regionInput: RIpOutput: ROp% Delete one inequality at a time% and check that region remains the sameROp = RIpforeach ineq ∈ RIp doR = ROp− region(ineq)if R == RIp thenROp = RendendFigure 4.6: Simplify Closed Region (x− 2y = 0 is removed).pass through any vertex.Example 4.4. In this example the blue region shown in Figure 4.7 is definedby four inequalities:x+ 7y + 4 ≤ 0,−x+ 95y + 5 ≤ 0,x+ 4 ≤ 0,y + 5 ≤ 0.The first inequality does not go through any feasible vertex and is eliminated.264.2. SOME METHODS FROM CLASS REGIONAlgorithm 4: Removing inequalities that do not pass through afeasible vertexInput: RIpOutput: ROpforeach inequality ∈ RIp doforeach vertex ∈ RIp doif inequality goes through vertex thenbreakendendmark inequality to be deletedendROp = RIp - region(marked inequalities)y + 5 x+7y+49y/5 - x + 5x+49y/5 - x + 5y + 5 x+4Figure 4.7: Remove inequalities that do not pass through a feasible vertex(x+ 7y + 4 = 0 is removed).274.2. SOME METHODS FROM CLASS REGIONSimplify unbounded region with linear and quadratic inequalities.This simplification process is done after applying Algorithm 4. Here wetry and remove redundant inequalities which go through valid vertices. Foreach vertex, we get all inequalities going through it. If only two inequalitiesgo through each vertex we already have a minimal set and nothing needs tobe done.For vertices with more than two polynomials going through them, we de-tect inequalities that can be removed. We consider two pairs of inequalities,the first pair would be the inequality being tested along with its adjacentinequality in the anticlockwise direction (all inequalities pass through thesame vertex so anticlockwise refers to the angle between the lines defined bythe inequality where ≤ is replaced by =) and the other pair would be the in-equality under consideration along with its adjacent inequality in clockwisedirection. For both pairs, we take the bisector line between the two linesinduced by the two inequalities. We can do this by plotting a line throughthe vertex with a slope between the slope of the two lines we are considering.This divides the region into four regions as shown in Figure 4.8. We test ifthese regions are feasible. A region is said to be feasible if every point inthat region satisfy the inequalities of the region we are trying to simplify sotesting any one point in the required region will be sufficient. If any regionis feasible, we cannot eliminate that inequality. If all the tested regions areinfeasible then the inequality can be removed. This process is explained inAlgorithm 5 and illustrated in the next example.Example 4.5. In Figure 4.9(a), we have four inequalities.y − 2 ≤ 0,x+ y − 6 ≤ 0,y ≤ 0,y2 − x ≤ 0.Three of these inequalities go through (4, 2). In Figure 4.9(b), the green linein the tangent to the parabola at the vertex (4, 2). The inequality y ≤ 2 isremoved after checking regions marked R1 . . . R8 to get the simplified regionin Figure 4.9(c).284.2. SOME METHODS FROM CLASS REGION9y/5 - x + 5y + 5 x+4R1R2R3R4(a) Before Simplificationx+49y/5 - x + 5(b) After SimplificationFigure 4.8: Process of simplifying an unbounded region with more than twoinequalities through a vertex.294.2. SOME METHODS FROM CLASS REGION(a) Before Simplification (b) Simplification Process(c) After SimplificationFigure 4.9: Simplifying a region with a quadratic inequality: Example 4.5.304.2. SOME METHODS FROM CLASS REGIONThis routine also detects singleton regions which can be removed.Example 4.6. Inequalities−1 + 52y + x ≤ 0,−x− 2y − 4 ≤ 0,x− 2y + 44 ≤ 0.all pass through (-24,10) but the region around it is infeasible. This is shownin Figure 4.10.314.2. SOME METHODS FROM CLASS REGIONAlgorithm 5: Simplifying an unbounded regionInput: RIpOutput: ROpforeach vertex ∈ RIp doStore the inequalities passing through it.endforeach inequality ∈ RIp doif inequality is linear thenStore the slope.endif inequality is quadratic thenStore the slope of the tangent.endend% slopes were computed to find adjacent inequalities.foreach vertex ∈ RIp doif number of inequalities through vertex = 2 thencontinueendforeach inequality through vertex doneighbor1 = inequality in anticlockwise direction.Get line between inequality and neighbor1.Take 4 points around vertex in the four regions createdbetween line, inequality and neighbor1.Test feasibility of the four regions.if feasible thencontinueendneighbor2 = inequality in clockwise direction.Repeat same procedure for neighbour 2Mark inequality for removal.endendRemove marked inqualities.324.2. SOME METHODS FROM CLASS REGION(-24,10)Figure 4.10: Singleton region.4.2.3 Set difference (minus)This operator is used to find the set difference of two feasible regions,R = R1 − R2. Output R is a list of regions. Algorithm 6 shows the steps;illustrative examples are given in figures 4.11 and 4.12.4.2.4 Maximum of two convex functions over a regionWe need to find the maximum of two functions over a given region inStep 1 and Step 3. In Step 1 the region is polyhedral and bounded. In Step3, the regions are parabolic and can be unbounded.When we find the maximum of two functions over a region, we eitherget one function over the given region or the region is partitioned into twopieces with one function maximum in one region and the other maximumover the second region.If there is a partition it is created due to the equality f1 = f2. Whilecomputing the maximum of two functions over a region, we need to know ifthere is a partition and if so which function is maximum over which region.Algorithm 7 computes and returns this information.In Step 1 the partition created is always polyhedral, while in Step 3 thepartition created is either polyhedral or parabolic. In Step 1 the region can334.2. SOME METHODS FROM CLASS REGIONAlgorithm 6: Computing the set difference between two regionsInput: R1, R2Output: Rif R1 == R2 thenR = ∅returnendl = list of inequalities from R1.v1 = end points of inequalities from R1.v2 = end points of inequalities from R2.rm = same inequalities in R1 and R2 with same end points.l2 = R2− rm. % Set difference.l = l − rm. % Set difference.l = l + (−l2). % Add inequalities (l2 ≥ 0) to the list l.R = region(l).if feasible (R) thenreturn.endCreate regions Ri with subsets of l.Store feasible Ri in a list R.return unique(R).344.2. SOME METHODS FROM CLASS REGION(a) R1. (b) R2. (c) R.Figure 4.11: Minus with one output region, R = R1 −R2be partitioned into two pieces only. In Step 3, the parabolic curve mightcreate more than two regions.As mentioned in the definition of Class region, we are assuming that theoriginal region has at least one vertex. Now if a partition is created, wewill get more vertices. This is crucial in computing the maximum using thisalgorithm,Algorithm 7 is called as[ l , ineq1 , ind ] = r.max(f1, f2)Here “l” is a logical flag which will be true if the region is not divided. Inthis case we return the function index of the maximum function in “ind” (1for f1 or 2 for f2). If there is a division, we return the inequality f1−f2 ≤ 0which is used to create one region and the negative of this inequality is usedto create the second region. The function which is maximum in the firstregion is returned in “ind”.354.2. SOME METHODS FROM CLASS REGIONAlgorithm 7: Computing the maximum of two functions over aregionInput: f1, f2, ROutput: l, ineq, ind% l returns true if one function is maximum over the entire region.% ineq returns the splitting inequality when l = false.% ind returns 1 or 2 depending on whether f1 or f2 is maximuml = trueif f1 == f2 thenind = 1returnendineq = f1 − f2R0 = R+ region(ineq)if R0 == ∅ ∥ R0 == R then% This means that the original region is not split% In this case we find which function is maximum as only% f1 or f2 can be maximum over the entire region.Choose point, P ∈ Rif f1(P ) > f2(P ) thenind = 1returnelseind = 2returnendend% ineq splits the regionl = falseChoose point, P ∈ R0if f1(P ) > f2(P ) thenind = 1returnelseind = 2returnend364.2. SOME METHODS FROM CLASS REGION(a) R1. (b) R2.(c) R3. (d) R4.Figure 4.12: Minus with two output regions (R1 −R2 = R3 +R4).37Chapter 5Convex envelope of eachpieceWe are given a PLQ function f . To compute the closed convex envelopeof a piece we utilize the method proposed in [Loc16]. We obtain a polyhedraldivision of the domain.Theorem 5.1 ([Loc16, Theorem 1.1]). For the bilinear and the fractionalterms of a function, the convex envelope over any polytope is characterizedby a unique polyhedral subdivision Pi : i ∈ I (possibly adding vertices; seeFigure 5.1). Over each member Pi, i ∈ I, the functional form of the convexenvelope is one of the following three forms:α6x2 + α5y2 + α4xy + α3x+ α2y + α1β3x+ β2y + β1, (5.1)α6x2 + α5y2 + α4xy + α3x+ α2y + α1, (5.2)α3x+ α2y + α1. (5.3)We obtain the convex envelope of a piece by solving the optimizationproblem given below. Given (x′, y′) ∈ P, the value of the convex envelope off over P is the optimal value of the optimization problem in three variablesa, b and c with an infinite number of constraints. The problem is formulatedas follows:conv(f + IP )(x, y) = maxa,b,cc, (5.4)subject tof(x′, y′)− [a(x′ − x) + b(y′ − y) + c] ≥ 0, ∀(x′, y′) ∈ P. (5.5)The optimal solution (a∗, b∗, c∗) of this problem gives a supporting hyper-plane for the convex envelope. In his paper [Loc16], Locatelli shows thatthe problem can be rewritten asconv(f + IP )(x, y) = maxa,b,cc, (5.6)38(a) Division with no new vertices. (b) Division with one added vertex.Figure 5.1: Polyhedral Subdivision.subject tof(xvi , yvi)− [a(xvi − x) + b(yvi − y) + c] ≥ 0, ∀i ∈ VP ,min(x′,y′)∈[vj1 ,vj2 ]f(x′, y′)− [a(x′ − x) + b(y′ − y) + c] ≥ 0, ∀j ∈ EP .where EP is the set of indices of the edges along which f is strictly convexand VP is the set of vertices which do not belong to any edge in EP . Foreach j ∈ EP , we have an edge Ej defined by its end points (vertices vj1 andvj2).We denote each edge by a line y′ = mjx′ + qj (as shown in Figure 5.2)and if we restrict f to this edge, we get f(x′) = f(x′,mjx′ + qj). Notethat if we have a vertical edge, the restriction of f to this edge gives us alinear function in y, which is not strictly convex; hence assuming the aboveequation for all such edges is not restrictive. Each vertex is denoted by(xvi , yvi). We now consider both cases.Case j ∈ EP : The constraint related to edge j ∈ EP with Ej = [vj1 , vj2 ]can be writtenminx′∈[xvj1 ,xvj2 ]{f(x′)− [a(x′ − x) + b((mjx′ + qj)− y)]} ≥ c, (5.7)where xvj1 ,xvj2 are the x-coordinates of the vertices vj1 and vj2 respectivelywith (xvj1 ≤ xvj2 ).Let the unconstrained minimum point of x′ 7→ fj(x′) − (a + bmj)x′ besj(a, b). We define sj(a, b) = ±∞ if the function is linear. If sj(a, b) ̸= ±∞,then f ′j(sj(a, b)) = a+ bmj .39Figure 5.2: Function restricted to edge.Let us now introduce three regionsD−j = {(a, b) : sj(a, b) ≤ xvj1},Dj = {(a, b) : xvj1 ≤ sj(a, b) ≤ xvj2},D+j = {(a, b) : sj(a, b) ≥ xvj2}.As f when restricted to the edge and is strictly convex, the first derivativeis increasing in the interval [xvj1 , xvj2 ], thus we haveD−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 )}.Therefore for each j ∈ Ep the optimal solution is (xj , yj), wherexj(a, b) =xvj1 if (a, b) ∈ D−j ,sj(a, b) if (a, b) ∈ Dj ,xvj2 if (a, b) ∈ D+j .and the associated y coordinate is yj(a, b) = mjxj(a, b) + qj . As xj(a, b) isa continuous function, we would get the same value on the common bound-aries.Case i ∈ VP : For each i in VP we have xi(a, b) = xvi and yi(a, b) = yvi .So if we define a function asηk(a, b) = f(xk(a, b), yk(a, b))− axk(a, b)− byk(a, b), ∀k ∈ VP ∪ EP , (5.8)then we can write equation (5.4) asconv(f + IP )(x, y) = max{c : ηk(a, b)+ax+ by ≥ c , ∀k ∈ VP ∪EP }. (5.9)40The convex envelope can then be computed by solving the convex prob-lem (5.9). Locatelli [Loc16] proved that at points different from vertices of Pand not belonging to edges where the function is strictly convex, there willbe at least two active constraints at an optimal solution. Using this fact, weproceed to solve the optimization problem as follow. We define subproblemsfor every pair (i, j), i, j ∈ EP ∪ VP , i ̸= j asmax{ηi(a, b) + ax+ by : ηi(a, b) = ηj(a, b),ηi(a, b) ≤ ηk(a, b), k ̸= i, k ̸= j,(a, b) ∈ Sr},(5.10)where Sr = ∩j∈EP τj , τj ∈ {D−j , Dj , D+j }, and j ∈ EP .Here r is the reindexing of the pair. Hence r = |EP |(i−1)+j, where (i, j)is a pair with both i and j ∈ Ep. This gives us a method to compute theconvex envelope. Locatelli [Loc16] further proved that the convex envelopeis characterized by a polyhedral subdivision of the polytope P .Note that the number of vertices of the polytope P may not be verylarge, all the subproblems have disjoint interiors of their feasible regions,and many of these subproblems can be discarded in advance when theyhave empty interiors.We split the subproblem (5.10) into many subproblems by adding con-straints (a, b) ∈ Sr, which allows us to uniquely identify the functional formof each function ηj for j ∈ EP . Here we will be adding the constraints(a, b) ∈ Sijr , which allows us to uniquely identify the functional form of ηiand ηj .Thus the subproblems we consider are as follows. If (i, j) ∈ EP , thenSijr = τi ∩ τj , where τi ∈ {D−i , Di, D+i } and τj ∈ {D−j , Dj , D+j }; the sub-problem we consider ismaxa,bηi(a, b) + ax+ by, (5.11)subject toηi(a, b) = ηj(a, b),ηi(a, b) ≤ ηj(a, b), k ̸= i, j,(a, b) ∈ Sijr .415.1. EXAMPLEFigure 5.3: Bilinear form with domain as defined in Example 5.1.5.1 ExampleWe now walk through the entire example of the bilinear form which isillustrated in Figure5.3.5.1.1 SubproblemsLet f(x, y) = xy defined over the polytope P = conv{V1, . . . , V4}, withvertices V1 = (0, 0), V2 = (2, 0), V3 = (2, 1) and V4 = (1, 1) as shown inFigure 5.3.Recall ηk(a, b) = f(xk(a, b), yk(a, b))−a(xk(a, b)−byk(a, b), ∀k ∈ VP∪EP .On the edge [V1, V4], y = x, f(x) = f(x, y) = xy = xx = x2, f ′(x) = 2x andf ′′(x) = 2. So on this edge, the function is strictly convex. This is the onlystrictly convex edge on the given Polytope P as the other edges are eitherx = c or y = c, so when the function is restricted to them, we get a linear425.1. EXAMPLEfunction.Therefore EP = {1}, E1 = [V1, V4] and VP = {2, 3}.Substituting values of xk(a, b), yk(a, b) in η, we get at V1, η(a, b) = 0 andat V4, η(a, b) = 1− a− b.For a convex edge, ηj(a, b) = −((a +mjb − qj)2)/(4mj) − bqj . For theedge [V1, V4], we have mj = 1, qj = 0. Thereforeηj(a, b) = −(a+ b)24.So associated with edge [V1, V4], we getη1(a, b) =0, a+ b ≤ 0,− (a+b)24 , 0 ≤ a+ b ≤ 2,1− a− b, a+ b ≥ 2.Similarly at V2 we get η2(a, b) = −2a and at V3, η3(a, b) = 2 − 2a − b.Now to get our subproblems, the possible pair of indices in VP ∪ EP are(1, 2), (1,3) and (2,3).For (1,2) and the region D+1 , we getmax 1− a− b+ ax+ by, (5.12)subject to1− a− b = −2a,1− a− b ≤ 2− 2a− b,a+ b ≥ 2.For (1,2) and the region D1, we getmax−2a+ ax+ by, (5.13)subject to−(a+ b)24= −2a,−2a ≤ 2− 2a− b,0 ≤ a+ b ≤ 2.Finally for (1,3) and the region D−1 , we get435.1. EXAMPLEmax ax+ by, (5.14)subject to0 = 2− 2a− b,0 ≤ −2a,a+ b ≤ 0.Steps i, ii, iiiThe subproblem (5.12) is of the type ηi, ηj linear. Eliminate variable b,using b = a+ 1, we getmax{a(x+ y − 2) + y : 12≤ a ≤ 1}, (5.15)and the optimal value for this subproblem isu1(x, y) ={12x+32y − 1 if x+ y ≤ 2,2y + x− 2 if x+ y ≥ 2.The subproblem (5.13) is of the type ηi(a, b) nonlinear, ηj(a, b) linearover Sijr . Setting z1 = a+ b and eliminating a using the first constraint, wegetmaxz218(−2 + x− y) + z1y, (5.16)z218− z1 + 2 ≥ 0,0 ≤ z1 ≤ 2.Since the quadratic inequality holds and −2+x−y ≤ 0 over P , the optimalvalue of this subproblem as a function of (x, y) isu2(x, y) ={12x+32y − 1 if x+ y ≥ 2,2y2y+2−x if x+ y ≤ 2.In subproblem (5.14), eliminating b, using b = 2− 2a, we getmax{a(−2y + x) + 2y : a ≤ 0, a ≥ 2}. (5.17)This problem has no feasible solutions.445.1. EXAMPLEFigure 5.4: Polyhedral Subdivision of domain for Example 5.1.Figure 5.5: Convex Envelope for Example 5.1.455.1. EXAMPLEStep ivThe convex envelope is given byconv(f + IP )(x, y) = max{u1(x, y), u2(x, y)}.Using u1 and u2 obtained in earlier step, we getconv(f + IP )(x, y) ={2y2y+2−x if x+ y ≤ 2,2y + x− 2 if x+ y ≥ 2,as shown in figure 5.5. The polyhedral subdivision of P isP1 = P ∩ {(x, y) : x+ y ≤ 2},P2 = P ∩ {(x, y) : x+ y ≥ 2},as shown in figure 5.4.Along with the convex envelope, we also store the individual componentsof the functional form as they are useful in finding the conjugate later. Thefunctional forms for the bilinear function are:r(x, y) =ξ1(x, y)2ξ2(x, y)+ ξ0(x, y), (5.18)r(x, y) =ξ1(x, y)2ζ+ ξ0(x, y), (5.19)r(x, y) = ζξ1(x, y) + ξ0(x, y), (5.20)where ξi(x, y) are linear functions in (x, y).The method to compute the convex envelope as given in [Loc16] is anexponential time algorithm but as our case applies to a specific domain wecan reduce the time complexity.Conjecture 5.2. We can find the convex envelope in linear time. In ourcase, our original domain is a simple polygon. So we can triangulate it,bounding the maximum number of convex edges as a function of the numberof vertices. It is possible to add vertices while solving the subproblems butexperimentally we add a maximum of two vertices. Hence we get a linearbound on the number of subproblems we solve.46Chapter 6Conjugate of each pieceIn this chapter we follow the method given in [KL20] to compute theconjugate of a convex bivariate rational function over a polytope. We obtainthe conjugate as a piecewise function as illustrated in Figure 6.1.After obtaining the convex envelope, we compute the conjugate of therational function over the polytope. We first compute the domain of theconjugate which turns out to be a parabolic subdivision. Then we find theconjugate expressions on that domain.It was shown in Chapter 5, that the closed convex envelope of a bivariatenonconvex quadratic function defined over a polytope is a piecewise ratio-nal function with a polyhedral subdivision. The conjugate of the convexpiecewise rational function, say r, is then computed asr∗ = maxj∈{1,2,...,nr}(rj + Idom rj )∗,where nr is the total number of pieces.6.1 Domain of the conjugateTo find the domain of the conjugate, we recall the following lemmasmentioned in [KL20].Lemma 6.1 ([BM23, Proposition 10.2]). For x ∈ Rn, a convex set P ⊂ Rn,and an indicator functionIP (x) ={0 if x ∈ P,+∞ otherwise.the subdifferential ∂IP (x) = NP (x) where NP (x) is the normal cone of P atpoint x.476.1. DOMAIN OF THE CONJUGATEFigure 6.1: Step 2: Computing the conjugate of a convex bivariate rationalfunction over a polytope.Lemma 6.2 ([Kum19, Lemma 4.2]). Consider a function g : Rn → R∪{∞},a polyhedral set P ⊂ Rn, and let f(x) = g(x) + IP (x), then∂f(x) = ∂g(x) +NP (x) =∂g(x) if x ∈ int(P ),∂g(x) +NP (x) if x ∈ bndry(P ),∅ otherwise.For PLQ functions, dom(f∗) = ∂f(dom(f)). The rational function ob-tained in the step above can be written asr(x, y) =ζ1(x, y)2ζ2(x, y)+ ζ0(x, y), (6.1)where ζi(x, y) are linear functions in (x, y).We now split the process of computing the domain of the conjugate bythe type of entities of the domain of the convex envelope.6.1.1 InteriorWe first compute the domain of the conjugate corresponding to the in-terior of the domain of the convex envelope.Proposition 6.3 ([Kum19, Proposition 4.3]). Let r be as defined in equation(6.1), then there exist a parabolic curve Cr(s) such that ∪x∈dom(r)∂r(x) ={s : Cr(s) = 0}.486.1. DOMAIN OF THE CONJUGATEFor the rational function r and a polytope P , define f(x) = r(x)+IP (x).Then for all x ∈ int(P ), the set ∪x∈int(P )∂f(x) is contained inside a parabolicarc.To compute the parabolic curve, we use Algorithm 8.Algorithm 8: Computing the parabolic curve from rational ex-pressionsInput: r expressionsOutput: curve% Algorithm that computes the convex envelope stores ξ1, ξ2 and ξ0% as symbolic expressions corresponding to% expressions given in (5.18), (5.19) and (5.20)% r(x) =ξ21(x)ξ2(x)+ ξ0(x), where% ξ0(x) = ξ01x1 + ξ02x2 + ξ00,% ξ1(x) = ξ11x1 + ξ12x2 + ξ10 and% ξ2(x) = ξ21x1 + ξ22x2 + ξ20% ∂q∂x1 =2ξ2(x)ξ1(x)ξ11−ξ21(x)ξ21ξ2(x)2+ ξ01% ∂q∂x2 =2ξ2(x)ξ1(x)ξ12−ξ21(x)ξ22ξ2(x)2+ ξ02% let t = ξ1(x)ξ2(x) then% ∂q∂x1 = 2ξ11t− ξ21t2 + ξ01% ∂q∂x2 = 2ξ12t− ξ22t2 + ξ02% dual variables (s1, s2) = (∂q∂x1, ∂q∂x2 )Extract ξ21, ξ22 and ξ20 (coefficients of ξ2(x))Extract ξ11, ξ12 and ξ10 (coefficients of ξ1(x))Extract ξ01, ξ02 and ξ00 (coefficients of ξ0(x))Create symbolic variable t%We have dual variables (s1, s2)vs1 = s1 − (2ξ11t− ξ21t2 + ξ01)vs2 = s2 − (2ξ12t− ξ22t2 + ξ02)% ξ22vs1 − ξ21vs2 cancels t2 term, hence we can solve for tvt = solve (ξ22vs1 − ξ21vs2, t )% substituting vt in vs1 gives us the parabolic curvecurve = subs(vs1,t, vt)496.1. DOMAIN OF THE CONJUGATE6.1.2 VerticesWe state the Lemma from [Kum19], in order to give the method we areusing to compute the domain of the conjugate corresponding to vertices.Lemma 6.4 ([Kum19, Lemma 4.10]). For r ∈ C1, P a polytope and v avertex. Let f(x) = r(x) + IP (x), then ∂f(v) is an unbounded polyhedral set.There is one vertex at which both numerator and denominator equalzero; the rational function is extended by continuity at that vertex. In orderto compute the domain of the conjugate in this case we use the followingconjecture.Conjecture 6.5 ([Kum19, Conjecture 4.11]). Let r be as in equation (6.1),f(x) = r(x) + IP (x) and v be a vertex of P with ζ1(v) = ζ2(v) = 0 . Then∂f(v) is a parabolic region.From Lemma 6.2, the subdifferential at the boundary is the sum of ∂g(x)and the normal cone at a point x. Hence for our computation, we begin bycomputing the normal cones for the vertices using Algorithm 9. Next at ver-tices where limits exist, we use Algorithm 10 to compute the subdifferential,∂f(x). Then we update the subdifferential at the vertices corresponding toConjecture 6.5 by continuity (Algorithm 11).Algorithm 9: Computing the normal Cone at each vertexInput: rOutput: NC% For each vertex we have two inequalities.foreach vertex Vi do% Compute the slope m of the normal% to the edge (Vi, Vi+1)% (x, y) are symbolic variables.NC(i,1) = y − Viy −m(x− Vix)% Set direction of inequality% s.t Vi+1 does not satisfy the inequalityif NC(i, 1)(Vi+1) < 0 thenNC(i,1) = - NC(i,1)end% Similarly set NC(i,2) as the normal% to the edge (Vi−1, Vi)end506.1. DOMAIN OF THE CONJUGATEAlgorithm 10: Computing the subdifferential at vertices wherelimits existInput: r, NCOutput: subdifferential, unDV% unDV is a logical array marking the vertices where limits do notexist(x, y): extract primal variables from rdrdx = ddx(r)drdy = ddy (r)[xlimitexists, limitdx] = limitvars→vertices drdx[ylimitexists, limitdy] = limitvars→vertices drdyfor i = 1: ∥vertices∥ doif limit does not exist thenunDv(i) = truecontinueendunDv(i) = falses1 = drdx(vertexi)s2 = drdy(vertexi)for j=1:2 dosubd(i,j) = (s1, s2) + NC(i,j)endend516.1. DOMAIN OF THE CONJUGATEAlgorithm 11: Computing the subdifferential at vertices wherelimits do not exist (by continuity)Input: r, subdifferential, unDVOutput: subdifferential % unDV is a logical array marking thevertices where limits do not existfor j = 1: number of vertices doif not unDv(i) thencontinueendem = j-1;if em == 0 thenem = nvendep = j+1;if ep == nv+1 thenep = 1endsubdifferential(j,1) = -subdifferential(em,1);subdifferential(j,2) = -subdifferential(ep,2);subdifferential(j,3) = curveend526.1. DOMAIN OF THE CONJUGATE6.1.3 EdgesWe now proceed to compute the domain of the conjugate correspondingto the edges of the polyhedron.Lemma 6.6 ([Kum19, Lemma 4.12]). For r ∈ C1, P a polytope, and anedge E = {x : x2 = mx1 + c, xl ≤ x1 ≤ xu}, between vertices xl and xu. Letf(x) = r(x) + IP (x), then⋃x∈ri(E)∂f(x) =⋃x∈ri(E){s+∇r(x) : s1 +ms2 = 0, s2 ≥ 0}.Proposition 6.7 ([Kum19, Proposition 4.13]). For r as in Equation (6.1),a polytope P, and an edge E ={x : x2 = mx1 + c, xv− ≤ x1 ≤ xv+}, betweenvertices v− and v+, let f(x) = r(x) + IP (x), then ∪x∈P∂f(x) is either aparabolic region or a ray.We use Algorithm 12 to compute the normal cones for the edges. Foreach edge we get two inequalities, one at each vertex such that the edge isincluded in the normal cone. Using the normal cone data, we can computethe subdifferential for the edges using Algorithm 10. Algorithm 13 thenadds the parabolic curve as an inequality to the conjugate domains of edgeswho have vertices where the derivative does not exist.Algorithm 12: Computing the normal cone at points on edges.Input: rOutput: NCforeach edgei(V1, V2) do% Compute the slope m of the normal% to the edgei at V1% (x, y) are symbolic variables.NC(i,1) = y − V1y −m(x− V1x)% Set direction of inequality% s.t V2 does not satisfy the inequalityif NC(i, 1)(V1) > 0 thenNC(i,1) = - NC(i,1)end% Similarly set NC(i,2) as the normal% at vertex V2end536.1. DOMAIN OF THE CONJUGATEAlgorithm 13: Computing the subdifferential at points on edgeswith vertices of 00 typeInput: r, subdE, unDV, curveOutput: subdE(x, y): extract primal variables from rdrdx = ddx(r)drdy = ddy (r)% The curve is precomputed% Here we only get the direction of the inequalityfor j = 1:number of vertices do% unDv(j) = true, when derivative exists at vertex jif unDv(j) thencontinueendif unDv(j+1) thencontinueend% vx is x coordinate of vertex j% vy is y coordinate of vertex j% We check direction of parabolic inequality.dv11 = drdx(vxj , vxj+1)dv12 = drdy(vxj , vxj+1)dv21 = drdx(vyj , vyj+1)dv22 = drdy(vyj , vyj+1)if curve((dv11 + dv21)/(2), (dv12 + dv22)/(2)) < 0 thencurve = -curveendsubdE(j,3) = curve;end546.2. CONJUGATE EXPRESSIONSV1V2V3E12E13(a) dom r(x1, x2). (b) dom r(x1, x2)∗.Figure 6.2: Conjugate Domain for Example 6.9From the above lemmas and corollary, [KL20] gives the following theo-rem.Theorem 6.8. Assuming r is as in equation(6.1), P is a polytope, andf(x) = r(x) + IP (x), then ∪x∈P∂f(x) admits a parabolic subdivision.Example 6.9. Givenr(x1, x2) =36x21 + 21x1x2 + 36x22 − 81x1 + 24x2 − 252−12x1 + 9x2 + 75 ,defined over a polytope with vertices v1 = (−1, 1), v2 = (−3,−3), v3 =(−4,−3); we compute the domain of the conjugate, and illustrate it in Fig-ure 6.2.6.2 Conjugate ExpressionsTheorem 6.10 ([Kum19, Theorem 4.24]). For r as in Equation (6.1), apolytope P, and f(x) = r(x) + IP (x), the conjugate f∗(s) has a parabolicsubdivision such that its expression over each member of its subdivision hasone of the forms (6.2).556.3. CONJUGATE - ALGORITHMIC APPROACHgf (s1, s2) =ψ1(s1, s2)ζ00√ψ 12(s1, s2)+ ψ0(s1, s2),gq(s1, s2) = ζ11s21 + ζ12s1s2 + ζ22s22 + ζ10s1 + ζ01s2 + ζ00,gl(s1, s2) = ζ10s1 + ζ01s2 + ζ00,(6.2)where ψ1, ψ0 and ψ 12are linear functions in s and ζij ∈ R.Similar to the domain, computation of conjugate expressions can bedivided into computation for each entity of P . Thus we compute the ex-pressions for vertices and edges while the subdifferential for the interior iscontained in a parabolic arc and hence the computation of the conjugate isdeduced by continuity. Algorithm 14 computes the conjugate expressionsat the vertices while Algorithm 18 computes the expressions at the edges.Algorithms 15, 16 and 17 are called by Algorithm 18. Refer to the appendixof [Kum19] for the derivations of Algorithms 15, 16 and 17.Algorithm 14: Computing the conjugate expression at verticesInput: r, unVOutput: expr% (s1, s2) are dual variable% (vxj , vyj) is the jth vertex% expr = sTx− r(x)for j=1: number of vertices doif unV(j) thenexprj = vxjs1 + vyjs2 − limit(x,y)→(vxj ,vyj)r(x, y)elseexprj = vxjs1 + vyjs2 − r(vxj , vyj)endend6.3 Conjugate - Algorithmic approachWe finally combine all the above steps into Algorithm 19 to compute theconjugate of one piece.566.3. CONJUGATE - ALGORITHMIC APPROACHAlgorithm 15: Computing the quadratic expression at edges whenthe coefficient of y is 0)Input: ξ0, ξ1, ξ2, cOutput: expr% ξ0(x) = ξ01x1 + ξ02x2 + ξ00,% ξ1(x) = ξ11x1 + ξ12x2 + ξ10 and% ξ2(x) = ξ21x1 + ξ22x2 + ξ20Extract ξ21, ξ22 and ξ20Extract ξ11, ξ12 and ξ10Extract ξ01, ξ02 and ξ00c2 =ξ202ξ12c3 = −ξ02c2c7 = ξ12c2d2 = −ξ12ξ02c2 + ξ13 − cξ11d3 = −ξ202c2 − cξ01 + ξ00c5 =c27ξ23c6 =2c7d2ξ20− c3d4 =d22ξ20+ d3expr = (c2 − c5)s22 − cs1 + (c3 − c6)s2 − d4576.3. CONJUGATE - ALGORITHMIC APPROACHAlgorithm 16: Computing quadratic expressions at edges whereξ21 +mξ21 = 0Input: ξ0, ξ1, ξ2,m, qOutput: expr% ξ0(x) = ξ01x1 + ξ02x2 + ξ00,% ξ1(x) = ξ11x1 + ξ12x2 + ξ10 and% ξ2(x) = ξ21x1 + ξ22x2 + ξ20Extract ξ21, ξ22 and ξ20Extract ξ11, ξ12 and ξ10Extract ξ01, ξ02 and ξ00t0 =−ξ01−mξ022(ξ11+mξ12)t1 =12ξ11+mξ12t2 =m2(ξ11+mξ12)γ10 = t1ξ23+qξ22ξ11+mξ12γ01 = t2ξ23+qξ22ξ11+mξ12γ00 =t0(ξ23+qξ22)−ξ13−qξ12ξ11+mξ12ζ11 = − (ξ11γ10+mξ12γ10)2ξ23+qξ22+ γ10ζ12 = −2∗(ξ11γ01+mξ12γ01)(ξ11γ10+mξ12γ10)ξ23+qξ22+γ01 +mγ10ζ22 = − (ξ11γ01+mξ12γ01)2ξ23+qξ22+mγ01ζ10 =−2(ξ11γ01+mξ12γ10)(ξ13+ξ11γ00+ξ12)(q+mγ00)ξ23+qξ22−mξ02γ10 + γ00 − ψ01γ10ζ01 =−(2(ξ11γ01+mξ12γ01(ξ10+ξ11γ00+ξ12(q+mγ00)))ξ23+qξ22−mξ02γ01 − ξ01γ01 +mγ00 + qζ00 =−(ξ13+xi11γ00+ξ12(q+mγ00))2ξ23+qξ22− ξ00 − ξ01γ00 − ξ02(q +mγ00)expr = ζ11s21 + ζ12s1s2 + ζ22s22 + ζ10s1 + ζ01s2 + ζ00586.3. CONJUGATE - ALGORITHMIC APPROACHAlgorithm 17: Computing the fractional expression at edges withξ21 +mξ21 ̸= 0Input: ξ0, ξ1, ξ2,m, qOutput: expr% ξ0(x) = ξ01x1 + ξ02x2 + ξ00,% ξ1(x) = ξ11x1 + ξ12x2 + ξ10 and% ξ2(x) = ξ21x1 + ξ22x2 + ξ20Extract ξ21, ξ22 and ξ20Extract ξ11, ξ12 and ξ10Extract ξ01, ξ02 and ξ00ζ00 = (ξ21 +mξ22)2δ1 = −2(ξ10ξ21 − ξ11ξ20 +mξ10ξ22 −mξ12ξ20−qξ11ξ22 + qξ12ξ21)(ξ01ξ21 + ξ211+m(ξ212 ∗m+ ξ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)ψ1 = 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ψ 12= −(ξ21 +mξ22)s1 −m(ξ21 +mξ22)s2 + ξ01ξ21+ξ211 +m(mξ212 + ξ01ξ22 + ξ02ξ21 + 2ξ11ξ12 +mξ02ξ22)ψ0 = (−ξ20(ξ21 +mξ22)− qξ22(ξ21 +mξ22))s1+(qξ21(ξ21 +mξ22)−mξ20(ξ21 +mξ22))s2 + δ0expr = ψ1ζ00√ψ 12+ ψ0596.3. CONJUGATE - ALGORITHMIC APPROACHAlgorithm 18: Computing the conjugate expression at edgesInput: ξ0, ξ1, ξ2Output: expr% Coded direct formulae from [Kum19, Appendix C]for j=1: nv (number of edges) do% Edge: ax+ by + c = 0Extract coefficients a,b,c% numbering starts from nv+1 as first nv are% expressions corresponding to verticesif b==0 thenexpr(nv + j) = conjugateExprQuad0(ξ0, ξ1, ξ2, c)elsem = −abq = − cbif ξ21 +mξ21 == 0 thenexpr(nv + j) = conjExprQuad(ξ0, ξ1, ξ2,m, q)elseexpr(nv + j) = conjugateExprFract(ξ0, ξ1, ξ2,m, q)endendend606.3. CONJUGATE - ALGORITHMIC APPROACHAlgorithm 19: Computing the conjugate of rational expressionsInput: rOutput: s (PLQ function)type = function.type % corresponding to (5.18), (5.19), (5.20)declare symbolic dual variables s1 and s2if type == 1 thenGet ξ0, ξ1 and ξ2 stored with convex envelopeGet parabolic curveGet Normal Cone at VerticesGet subdifferential at VerticesUpdate subdifferential for vertices where limits do not existGet Normal Cone at EdgesGet subdifferential at EdgesUpdate subdifferential for edge where limits do not existGet conjugate expressions for verticesGet conjugate expressions for edgesend% type 2 not implemented as it is a subset of type 1if type == 3 thenGet ξ0Get Normal Cone at VerticesGet subdifferential at VerticesGet conjugate expressions for verticesendStore domain and corresponding expressions as si61Chapter 7Conjugate of a PLQ functionWe started with a bivariate PLQ function,f = mini=1...nq(qi + Idom qi).In Step 1, we found the convex envelope of each piece. This gave us apolyhedral subdivision of the domain of each piece,conv(qi + Idom qi) = minj=1...nri(ri,j + Idom ri,j ).See Figures 3.3(b) and 3.3(c) for an illustrative example.In Step 2 we then found the conjugate of each of these convex envelopes,(ri,j + Idom ri,j )∗ = mink=1,...nsi,j(si,j,k + Idom si,j,k),which have a parabolic subdivision with full domain; it is illustrated inFigures 3.3(d)–(f).In this chapter we explain Step 3 where we will find the conjugate f∗ ofthe original PLQ function f by computing the maximum of all the conjugatesobtained in Step 2.7.1 Method for Step 3 of the main algorithm.In order to find the maximum of the conjugate of the entire PLQ func-tion, we work in two stages as follow (see Figure 7.1).Step 3a: Find the maximum of the conjugate for each piece i. The maximumfor each piece i is a piecewise function denoted by si,k such thatmaxj=1...nrimink=1,...nsi,j(si,j,k + Idom si,j,k) = mink=1...nsi(si,k + Idom si,k).Step 3b: Find the maximum of the conjugate across these pieces. This maxi-mum is again a piecewise function that we denote as sk so thatmaxi=1...nq(si,k + Idom si,k) = mink=1...ns(sk + Idom sk) = q∗.627.1. METHOD FOR STEP 3 OF THE MAIN ALGORITHM.(a) Step 3a: Conjugate of one piece.(b) Step 3b: Conjugate of PLQ function.Figure 7.1: Step 3: Computing the conjugate of a PLQ function.637.2. STEP 3A7.2 Step 3aIn Step 3a, we find the maximum of the conjugate for one piece. Asstated earlier, on finding the convex envelope of a piece, we get rationalfunctions rij defined on a polyhedral subdivision. Corresponding to eachrij , we find a conjugate si,j,k, where ∪k dom si,j,k = R2. Thus for each piecewe have many piecewise conjugates. To find the conjugate for one piece weneed to find the maximum of all these conjugates.We divide the task of finding the maximum of these conjugates into threesteps. In the first step (i) we consider only the domains, in the second step(ii) the functions defined on the domains, and in the third step (iii) we mergeadjacent regions if they have the same conjugate expressions.7.2.1 Step 3a(i)Let us first consider the domains of the conjugates corresponding to thefirst two convex functions in the ith piece. The domain of the first conjugateis ∪k1=1...nsi,1 dom si,1,k1 while ∪k2=1...nsi,2 dom si,2,k2 is the domain of thesecond conjugate of the ith piece. We need to find the intersection betweenthese domains, dom sti,k = dom si,1,k1 ∩ dom si,2,k2 , for all k1, k2.Currently we are exhaustively finding the Cartesian product betweentwo sets of regions and checking if the intersection of the interior of thetwo regions is nonempty. We use the Algorithm 2 to find the intersection.This gives a complexity of O(n2), where n = max(nsi,1 , nsi,2), but mightbe done in O(n) as mentioned in conjecture 7.1 below. When we find anonempty intersection, we store it in the output list along with the twofunctions which were defined on the two regions we were intersecting. Sowe are basically storing nonempty intersections along with two functions,as (dom sti,k, [si,k,1, si,k,2]) = (dom si,1,k1 ∩ dom si,2,k2 , [si,1,k1 , si,2,k2 ]). Theoutput dom sti,k is a list of nonempty domains covering the entire R2 planeas shown in Figure 7.2.Algorithmic Approach: R = R1 ∩R2In Algorithm 20 we are finding the Cartesian product of two sets offeasible regions. This routine assumes that each set covers the entire planeR2. We traverse the two sets and take two regions, one from each andthen find their intersection as explained in Algorithm 2. We then simplifythe region obtained as explained in Algorithm 4. After simplification if theintersection has a nonempty interior, we store it in the output set.647.2. STEP 3A123456(a) dom si,1,k1 , k1 = 1 . . . 6.123(b) dom si,2,k2 , k2 = 1 . . . 3.(1,1) (2,1)(2,3)(3,1)(3,4)(4,3)(5,1)(5,2)(6,2)(6,3)(c) dom sti,k = dom si,1,k1 ∩ dom si,2,k2 .Figure 7.2: Intersection of domains of two conjugates.657.2. STEP 3AAlgorithm 20: Computing the Cartesian product of two lists ofregionsInput: R1, R2Output: Rn=0foreach (r1, f1) ∈ R1 doforeach (r2, f2) ∈ R2 dorf = r1 ∩ r2if rf == ∅ thencontinueendrf = rf.simplifyUnboundedRegion;if rf == ∅ thencontinueendn = n + 1R(n) = (rf, [f1, f2]);endendAlthough Algorithm 20 is computed in O(n2), we can reduce its com-plexity.Conjecture 7.1. Step 3a(i) can be computed in O(n) worst-case time, wheren is the sum of the entities of the regions considered.When we consider regions that our neighbours. We get nonempty inter-sections only when the original convex envelopes have a common vertex oredge. So a propagation algorithm could provide a linear time algorithm.Algorithm 21 computes the intersection in O(n). We have explained thisalgorithm with the help of Example 7.2 which is illustrated in Figure 7.3.Example 7.2. Figure 7.3(a), shows the convex polyhedral subdivision ofthe domain of one piece as obtained in Step 1. As a result of Step 1 we getri1 = conv(V1, V2, V3) and ri2 = conv(V2, V3, V4).Figure 7.3(b) and Figure 7.3(c) show dom(si1) and dom(si2) respec-tively, where sij = r∗ij . From Figure 7.3(a) we observe that dom(ri1) ∩dom(ri2) = {V2, V3, E23}. The set dom(si1) has pieces corresponding toR1 = {V1, V2, V3}, whereas the dom(si2) has pieces corresponding to R2 ={V2, V3, V4, E23, E24, E34}.667.2. STEP 3AWe can split the entities in R1 and R2 based on the intersection weobserved in the primal. Thus R1 can be split as RC1 ∪ RC ′1, where RC1 ={V2, V3} and RC ′1 = {V1} and R2 = RC2 ∪RC′2 ,where RC2 = {V2, V3, E23}and RC′2 = {V4, E24, E34}. Let SCj be the pieces in dom(sij) correspondingto RCj . Then we get dom(si1) ∩ dom(si2) ⊆ (SC1 ∩ SC2) ∪ (SC ′1 ∩ SC2) ∪(SC1 ∩ SC ′2).In Algorithm 21, we have two sets R1 and R2 which have the entitiesfrom the domains of the convex envelope. First we find R1 ∩ R2. This canbe done by traversing R1 and marking the entities in a logical array. Thenwe traverse R2 and check if any entities are marked and put them in theintersection if marked. This intersection takes number of operations equalto the sum of entities in R1 and R2. Then we have three Cartesian products.The intersectionR1∩R2 has at most three entities (two vertices and an edge),thus the first Cartesian product takes 9 operations. The second and thirdtake three times the number of remaining entities resulting in maximumnumber of operations of 4 · (sum of entities).Algorithm 21: Computing the intersection of two lists of dualregions that are adjacent in the primal domainInput: R1, R2, S1, S2Output: SWhen conjugates are computed - store entity along with each piece.R = R1 ∩R2.SC1 = S1(R)SC2 = S2(R)SC ′1 = S1 − SC1SC ′2 = S2 − SC2S = (SC1 SC2) ∪ (SC ′1 SC2) ∪ (SC1 SC ′2)DivisionsAt this stage we are finding the intersection of domains of the originalconjugates. These domains are parabolic subdivisions [KL20]. After com-puting the intersection we again get a parabolic subdivision.Proposition 7.3. The output of Step 3a(i) is a parabolic subdivision.Proof. Each dom si,j,k is a parabolic subdivision, which is the nonemptyinterior of a set of parabolic inequalities. So each domain is of the form677.2. STEP 3A(a) conv(qi + Idom qi) =minj=1...2(rij + Idom rij ).(b)dom(si,1,k).V1V3V2(c)dom(si,2,k).E34V4E24E23V2V3(d) dom(si,1,k) ∩ dom(si,2,k).Figure 7.3: Computing dom(si,1,k)∩dom(si,2,k) using dom(ri,1)∩dom(ri,2).687.2. STEP 3Adom si,j,k = {(x, y), alx2+ blxy+ cly2+ dlx+ ely+ fl < 0, b2l − 4alcl = 0, l =1 . . .K} whereK ∈ N is the number of inequalities defining dom si,j,k. Whenwe find the intersection of two such domains, we only store the regions withnonempty interiors. These are again defined by parabolic inequalities of thesame form, hence the subdivision is still parabolic.7.2.2 Step 3a(ii)In Step 3a(i) we get a list dom sti,k, along with two functions si,k,1 andsi,k,2 defined over each dom sti,k. We now find the maximum function sti,k =max{si,k,1, si,k,2} over the given dom sti,k. In doing so it is possible thatdom sti,k is further divided into two regions one where si,k,1 is the maximum,and the other where the maximum is si,k,2; or we get the same region witheither si,k,1 or si,k,2 as the maximum.Example 7.4. The blue region in Figure 7.4 is denoted by B while R is thered region.B = {(u, v) : −u− 7v − 4 ≤ 0,u+ 7v − 10 ≤ 0,−u− 2v − 4 ≤ 0,u+ 2v − 4 ≤ 0,48u− 56v + 4uv + u2 + 4v2 − 184 ≤ 0}.R = {(u, v) : −u− 7v − 4 ≤ 0,148u− 196v + (u+ 7v)2 − 684 ≤ 0,u+ 2v − 4 ≤ 0,56v − 48u− 4uv − u2 − 4v2 + 184 ≤ 0}.Functions f1(u, v) = −5u+5v+25 and f2(u, v) = u28 + uv2 +u+ v22 −2v+2are defined on B ∪ R. When we find the max{f1, f2} over B ∪ R, we geta further subdivision, given by the equality f1 − f2 = 0. This gives us thesubdivision, B and R as shown in Figure 7.4. Function f1 is the maximumover B and f2 over R.697.2. STEP 3AFigure 7.4: Further division by a parabolic curve as explained in Example7.4.Algorithmic Approach: Maximum of two functions over a givendomain.This algorithm gives us a way to compute the maximum of two functionsover a region and split the region if required. The merge routine, definedas Step 3a(iii), is called at the end of this function. The routine to find if afunction is maximum over a region is given in Algorithm 7.707.2. STEP 3AAlgorithm 22: Computing the maximum of two functions over agiven domainInput: R1Output: R2foreach (r, [f1, f2]) ∈ R1 do[l, ineq, i] = maximum(r, [f1, f2])% l = true when fi is maximum over rif l thenn = n + 1R2(n) = (r, find)continueendr1 = r + region(ineq)r2 = r + region(−ineq)r1 = simplifyUnboundedRegion(r1) %Algorithm 5r2 = simplifyUnboundedRegion(r2)if ind == 1 thenn = n + 1R2(n) = (r1, f1)n = n + 1R2(n) = (r2, f2)elsen = n + 1R2(n) = (r1, f2)n = n + 1R2(n) = (r2, f1)endendmerge717.2. STEP 3ADivisionsAt this point, when we find the maximum of the two functions definedon the intersecting domain, we get a further subdivision with an inequalityf1−f2 ≤ 0 for one region and f2−f1 ≤ 0 for the other. Thus in the maximumwe can get a subdivision which is the difference of conjugate functions. From[KL20, Theorem 4.24, 4.27] we get the conjugate functions of the forms (6.2).In order to find the possible subdivisions, we enumerate all the pairs, of theconjugate functions, (f1, f2) that give us a further subdivision of the domain.We first give propositions for each possible case and the partition it can giveus and then put everything together in Theorem 7.9.Proposition 7.5. Let the two functions defined on the parabolic domain belinear. If we find the max{f1, f2}, we can get a further subdivision of thedomain. This division is given by a linear inequality hence the subdivisionremains parabolic.Proof. The domain is split by the equality, f1 − f2 = 0 which is linear.In order to get the subdivision when one conjugate expression is quadratic,we first need to show that the quadratic expression gq obtained in (6.2) isalways parabolic.Proposition 7.6. The quadratic expression gq obtained in (6.2) is parabolic.Proof. The conjugate expression obtained is a quadratic only in a certaincondition. In Step 1, when we find the convex envelope of each piece, wecan obtain a bivariate rational function,ξ21(x)ξ2(x)+ξ0(x), where ξ0(x) = ξ01x1+ξ02x2+ξ00, ξ1(x) = ξ11x1+ξ12x2+ξ10 and ξ2(x) = ξ21x1+ξ22x2+ξ20. Thenin Step 2, we find the conjugate expression of this rational function. Usingthe method given in the paper [Kum19], when we find the expression on anedge (y = mx + q) of the polytope and ξ21 +mξ22 = 0, we get a quadraticexpression. The quadratic expression obtained isf∗(s) = ζ11s21 + ζ12s1s2 + ζ22s22 + ζ10s1 + ζ01s2 + ζ00, (7.1)where,ζ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,727.2. STEP 3Aζ22 = −(ξ11γ10 +mξ12γ10)2ξ20 + qξ22+ γ10m,ζ10 = −2(ξ11γ10 +mξ12γ10)(ξ10 + ξ11γ00 + ξ12(q +mξ00))ξ20 + qξ22+γ00 −mξ02γ10 − ξ01γ10,(7.2)ζ01 = −2(ξ11γ01 +mξ12γ01)(ξ10 + ξ11γ00 + ξ12(q +mξ00))ξ20 + qξ22+mγ00 −mξ02γ01 − ξ01γ01 + q,(7.3)andζ00 = −(ξ10 + ξ11γ00 + ξ12(q +mξ00))2ξ20 + qξ22− ξ00 − ξ01γ00 − ξ02(mγ00 + q).Now substituting the above values, we obtainζ212 − 4ζ11ζ22 = 0.Hence (7.1) is parabolic.The MATLAB code to verify this proof is given in Appendix B.Lemma 7.7. If f1 is a parabolic quadratic function and f2 is linear, com-puting max{f1, f2} may result in a division by the curve f1−f2 = 0. In thatcase, the curve f1 − f2 = 0 is parabolic.Proof. Let the parabolic expression defined be f1(s) = ζ111s21 + ζ112s1s2 +ζ122s22+ζ110s1+ζ101s2+ζ100 and the linear function be f2(s) = ζ210s1+ζ201s2+ζ200.The new subdivision is given by f1(s) − f2(s) = ζ111s21 + ζ112s1s2 + ζ122s22 +(ζ110− ζ210)s1+(ζ101− ζ201)s2+(ζ100− ζ200) = 0. As f1(s) is parabolic, we haveζ212 − 4ζ11ζ22 = 0, hence f1(s)− f2(s) is also parabolic.This is shown in Figure 7.4.Proposition 7.8. When one function is fractional and the other is linear,the domain will be divided by a parabolic curve.Proof. Notef1(s) =ψ1(s1, s2)ζ00√ψ2(s1, s2))+ ψ0(s1, s2), (7.4)where ψ0(x) = ψ01x1 + ψ02x2 + ψ00, ψ1(x) = ψ11x1 + ψ12x2 + ψ10 andψ2(x) = ψ21x1 + ψ22x2 + ψ20; andf2(x) = ψ31x1 + ψ32x2 + ψ30. (7.5)737.2. STEP 3AThe fractional expression is given by Expression (7.4), while the linearexpression is given by Equation (7.5). Now the subdivision is given byf1(s)− f2(s) = 0, which can be simplified as follows:ψ1(s1, s2)ζ00√ψ2(s1, s2))+ ψ0(s1, s2) = f2(s1, s2),ψ1(s1, s2)ζ00√ψ2(s1, s2))= f2(s1, s2)− ψ0(s1, s2),ψ1(s1, s2)2ζ200(ψ2(s1, s2))= (f2(s1, s2)− ψ0(s1, s2))2,ψ1(s1, s2)2 = ζ200(f2(s1, s2)− ψ0(s1, s2))2ψ2(s1, s2).This is a degree 3 equality. But when you substitute symbolic values ofψ2(s1, s2), ψ1(s1, s2), ψ0(s1, s2) and ζ00 from Chapter 6, the degree 3 termsvanish. Then by taking coefficients of s21, s22 and s1s2 and substitutingsymbolic values you get ζ212−4ζ11ζ22 = 0. Hence the new inequality is againparabolic.We have written a code in MATLAB to verify the above proof by sub-stituting symbolic expressions and simplifying; see Appendix B.Theorem 7.9. Computing the maximum of two conjugates gives us a parabolicsubdivision.Table 7.1: Possible divisions of domain obtained in Step 3a (ii).si,1,k si,2,k si,1,k − si,2,kLinear Linear LinearLinear Parabolic ParabolicLinear Fractional ParabolicProof. When we find the convex envelope of each piece in Step 1, we donot get rational functions in adjacent domains. Rational functions ariseonly in (D.7). In this case, the region is divided into three regions. Themiddle region has the rational function, while the adjacent regions havelinear expressions. In Step 2, only the rational functions restricted to an747.2. STEP 3AFigure 7.5: Merging adjacent pieces with the same function.edge give rise to quadratic or fractional expressions. In Step 3, when wecompare two functions we always get one of them as linear. Nonemptyintersections are obtained only when the primal has a common edge. Hencethe two functions we compare always have one linear function. Thus theadditional subdivision is one of the forms given in table 7.1. The proofsthat these subdivisions are parabolic are given in propositions 7.5, 7.7 and7.8.7.2.3 Step 3a(iii)At this point we computed the maxima and obtained the conjugate asa piecewise function. This can be further simplified by merging adjacentregions where functions are the same as shown in Figure 7.5. Coefficients ofthe functions involved are stored as rational numbers, so we can determinetwo identical functions without floating point considerations.The merge routine is called at the end of Algorithm 22.7.2.4 Overall algorithm for Step 3a.For every piece of the PLQ function, we have multiple conjugates, rijk.So far, we have given a method to compute the maximum of two conjugates.We can now add one conjugate at a time to get the conjugate of the entirepiece as given in Algorithm 23. We first copy the first conjugate into the757.3. STEP 3Boutput. We then iterate over the remaining pieces, performing steps 3a(i)and 3a(ii) and finally get the conjugate of the entire piece.Algorithm 23: Computing the maximum to obtain the conjugateof a pieceInput: POutput: Pmax%Copy the first conjugatePmax(d, f) = P1.conjugate(d, f)foreach i ∈ (2 . . . |P |) do% Find intersection of domainsPmax(d, f1, f2) = Pmax(d, f) Pi.conjugate(d, f2)% Find maximum function on each domainPmax(d, f) = Pmax(d, f1, f2).maximumPendDivisionConjecture 7.10. The function maxj(rijk) admits a parabolic subdivision.By Proposition 7.3 and Theorem 7.9, the maximum of two conjugatesgives us a parabolic subdivision. In computing the convex envelope of a pieceonly one division can have a rational function defined on it. In computing theconjugate, only the rational function would give us a quadratic or fractionalform. Thus computing the maximum of multiple divisions would give us aparabolic subdivision.Example 7.11. In the example in Figure 3.3(a), we had two pieces. Onfinding the convex envelope in Step 1 we obtained convex functions definedover a polyhedral subdivision, refer to Figure 3.3(b-c). Corresponding tothese, in Step 2 we got two piecewise conjugates as shown in Figure 3.3(d-g). The maximum of the conjugates for Piece 1, is shown in Figure 7.6.The same process is repeated for Piece 2, shown in Figure 3.2 and themaximum conjugate for this piece is shown in Figure 7.7.7.3 Step 3bIn Step 3a, corresponding to each piece of the PLQ function, we obtaineda maximum conjugate which is the maximum of all the conjugates for that767.3. STEP 3B(a) s1,1,k. (b) s1,2,k.(c) max{s1,1,ks1,2,k}.Figure 7.6: Maximum Conjugate for piece 1.777.3. STEP 3B(a) s2,1,k. (b) s2,2,k.(c) max{s2,1,ks2,2,k}.Figure 7.7: Maximum Conjugate for piece 2.787.3. STEP 3Bpiece. In this stage we find the maximum over all these maximums to getone maximum for the entire PLQ function. The steps followed to computethis are very similar to those followed in Step 3a. In Step 3b(i) we findthe intersection of the domains, dom stk = dom s1,k ∩ dom s2,k and in Step3b(ii), we find the maximum of two functions over the same domain asstk = max{sk,1, sk,2} over dom stk and then merge to get the maximum overthe first two pieces. We then iterate Step 3b(i) and Step 3b(ii), over theremaining pieces, dom stk = dom stk∩dom sj,k, stk = max{stk, sj,k} and mergeto finally get the maximum conjugate sk = stk, which is the conjugate ofthe PLQ function. This process is given in Algorithm 24. Algorithm 24 issimilar to Algorithm 23 but it is working on the conjugates of entire pieces.Algorithm 24: Computing the maximum to obtain the conjugateof the PLQ functionInput: PLQInOutput: PLQOut%Copy max conjugate of first piece% PLQIni refers to the ith piecePLQOut(d, f) = PLQIn1.maxC(d, f)foreach i ∈ 2 . . . |PLQIn| doPLQOut(d, f1, f2) = PLQOut(d, f) PLQIni.maxC(d, f2)PLQOut(d, f) = PLQOut(d, f1, f2).maximumPend7.3.1 DivisionsIn this step we are finding the maximum between two pieces. In Step3b(i), we find the intersection of parabolic subdivisions and hence get parabo-lic subdivisions as given in Proposition 7.3. In Step 3b(ii), we are finding themaximum of two functions over a parabolic domain. The pieces are adjacentto each other in the primal and we conjecture that we would get a parabolicsubdivision as if we had considered both pieces together as one piece.Example 7.12. Again, consider the example given in Figure 3.2. We hadtwo pieces. Corresponding to each piece, in Step 3a we got the maximumconjugate of each piece which were shown in Figures 7.6 and 7.7. We nowfind the maximum over these maximums to get the conjugate of the PLQfunction as shown in Figure 7.8.797.3. STEP 3B(a) s1,k. (b) s2,k.(c) max{s1,ks2,k}.Figure 7.8: Maximum Conjugate over all pieces.807.3. STEP 3BExample 7.13. Here we refer to the PLQ function shown in Figure 3.2 butinstead of dividing it into two pieces we take just one piece as illustrated inFigure 7.9 and use this to test that the final output for both examples is thesame. In this example as there is only one piece, Step 3b is not required.Observing the conjugates we computed in figures 7.9, 8.1 and 8.4, wecan estimate the number of pieces in the conjugate.Conjecture 7.14. The number of pieces in the conjugate (computed as amaximum) is equal to the number of vertices plus the number of convex edgesof the overall domain.817.3. STEP 3B(a) q1.(b) r1,j .(c) s1,1,k. (d) s1,2,k. (e) s1,3,k.(f) sk.Figure 7.9: Conjugate of the PLQ function from Example 7.13.82Chapter 8ExperimentsWe performed a lot of experiments to test and time the code. The timingsare recorded on a Intel(R) Xeon(R) CPU E5-1620 v4 @ 3.50GHz 3.50 GHzmachine with 64.0 GB RAM, We have run the code with MATLAB R2023Bon a Windows 10 operating system. The code is not optimized and wouldbe much faster in practice. Some of the results are given in this chapter.8.1 VerificationIn order to verify the code works correctly, we first verified small prob-lems by solving them by hand. Next we computed conjugates for examplesfrom [Kum19] and verified that we obtain the same results. Finally we di-vided the same domain in different ways and check that we obtained thesame conjugate (in the first two steps, the convex envelopes and conjugatescomputed are different). In all cases, we found that Step 3 computation iscorrectly implemented.8.1.1 Experiment 1We start with a function = xy defined on a domain = conv{V1, V2, V3, V4}where V1 = (0, 0), V2 = (1, 1), V3 = (2, 1) and V4 = (2, 0). This is theexample given in [Kum19].This function and the domain of its conjugate are shown in Figure 8.1.We repeat this experiment twice, once with one piece, second time with twopieces and obtain the same conjugate.8.1.2 Experiment 2We divide the same domain given in Example 3.3 in three different ways:1. The first example has two pieces and is illustrated in 8.2(a).2. The second example has one pieces and is illustrated in 8.2(b).838.1. VERIFICATION(a) Domain 1 (b) Domain 2 (c) Function(d) Domain of the conjugateFigure 8.1: Experiment 1: Computing the conjugate of a function with thedomain stored as one piece gives the same result as the same function withthe domain stored as two pieces.848.2. TIMING TESTS3. The third example again has two pieces and is illustrated in 8.2(c).Polyhedral division for this example is given in 7.8(c).Although the division is different, all three examples give the same con-jugate as illustrated in 7.13.8.2 Timing TestsWe find the conjugate of xy over a square conv{V1, V2, V3, V4} whereV1 = (−1,−1), V2 = (−1, 1), V3 = (1, 1) and V4 = (−1, 1). We then dividethis piece into n2 pieces and repeat the experiment where n = 1, 2, 3, 4 andtabulate the results. All three tests give the same output which is given inFigure 8.4(a) and the graph of the time to find conjugate vs. number ofpieces and conjugate vs edges is plotted in Figure 8.4(b).The results for n = 1, n = 2 and n = 3 are tabulated in Tables 8.1,8.2 and 8.3. The tables contain the computation time for each step foreach piece, and the average time per piece, which allow to estimate theimprovement we could get by parallelisation of the code. The results aresummarized in Table 8.4. We plot graphs for time against the number ofpieces originally given and the number of edges of the original domain. Boththese graphs appear quadratic.Table 8.1: Computational time in seconds: Number of Pieces = 1.convex pieces 2Time for Step 1 60.3Time for Step 2 3.6Time for Step 3a 9.1Time for Step 3b 0.0Total 73.1Using a better algorithm, we can reduce the computation time to linear.Conjecture 8.1. Computing conjugate of a PLQ function is linear withrespect to number of edges.Timings for both Step 1 and Step 3 can be improved significantly. InStep 1, a lot of subproblems can be discarded in advance. This code has notbeen implemented. In Step 3, we compute the Cartesian product of two listwhile computing the intersection of domains. The time taken for these loopscan be reduced as this code can be implemented in linear time. Merging858.2. TIMING TESTS(a) Domain 1 (b) Domain 2(c) Domain 3Figure 8.2: Domains used in Experiment 2.868.2. TIMING TESTS(a) r1. (b) r2.Figure 8.3: Polyhedral subdivision of domain of convex envelope for 8.2(c).(a) Domain of the conjugate in allcases.Time Vs PiecesTime Vs Edges(b) Time to compute the conjugatevs. pieces and vs. edges.Figure 8.4: Conjugate of xy over unit square with each piece subdividedinto n2 pieces and corresponding computation time per piece and per edge.878.2. TIMING TESTSTable 8.2: Computational time in seconds. Number of Pieces = 4; eachpiece is split in 2 convex pieces.Piece No 1 2 3 4 total averageStep 1 57.9 66.5 66.2 55.2 245.8 61.5Step 2 5.3 4.9 4.8 2.9 17.9 4.5Step 3a 12.6 12.7 12.9 8.0 46.2 11.6Average per piece 77.5Step 3b 168.4Total 478.3Table 8.3: Computational time in seconds. Number of Pieces = 9; eachpiece is split in 2 convex pieces.Piece No 1 2 3 4 5 6 7 8 9 total avgStep 1 48.4 46.8 48.6 49.6 47.8 46.9 46.9 47.7 49.1 431.8 47.9Step 2 3.4 3.7 3.6 3.5 3.6 3.9 3.7 3.5 3.6 32.4 3.6Step 3a 9.1 9.9 3.6 9.1 8.9 11.2 9.2 9.1 9.6 79.5 8.8Average 60.4Step 3b 1,661.4Total 2,205.1Table 8.4: Timing summarynPieces nVertices nEdges avg time for 1 piece total time1 4 4 73 1464 9 12 77 4789 16 32 60 2,20516 25 40 58 6,803888.2. TIMING TESTSadjacent regions where possible would also decrease list sizes, thus giving usa faster time.89Chapter 9ConclusionIn this thesis we have given a method to find the conjugate of a bivariatepiecewise linear-quadratic function. We have implemented the first threesteps of this method to get the conjugate of a bivariate PLQ function. Itis possible to extend the code to implement Step 4 to find the biconjugatein order to get the convex envelope of the PLQ function. In the currentimplementation, the complexity of Step 1 can be exponential in the worst-case, Step 2 is linear and Step 3 is polynomial. But we expect the complexityto be reduced in the future.9.1 Future WorkThe current implementation is in MATLAB. The code is written in away that could be ported to other languages and architectures.9.1.1 BiconjugateNow that we have the conjugate of the PLQ function, the next naturalstep is to compute the biconjugate in order to obtain the convex envelope.To compute the biconjugate, we need to compute the conjugate of frac-tional forms on parabolic subdivisions.9.1.2 Step 1As we have implemented this computation in four steps, we can improveeach of these steps. In Step 1, we can reduce the complexity of the algorithm.When there is no convex edge - complexity is(n2). When there is a convexedge, with a vertex opposite it - this gives a subdivision which is generallythe maximum on solving the subproblems. Thus when there is a convexedge, we can reduce the number of subproblems being solved.909.1. FUTURE WORK9.1.3 Step 3Step 3a can also be improved and implemented with a lower complexity.We have given an algorithm to implement the intersection of two domainsin O(n). It has also been observed that the conjugates corresponding tovertices and edges that are not a part of the division when we find the convexenvelope of each piece appear exactly the same in the final conjugate. Thisfact can be used to develop a better algorithm to find the maximum. If youobserve Figure 9.1(a), the conjugates corresponding to the three vertices(V2, V3, V6) as illustrated in Figure 9.1(b-c), which are not a part of thesubdivision, appear exactly the same in the maximum 9.1(d).9.1.4 Parallel ComputingWe could develop a parallel algorithm to find the convex envelope ofbivariate functions over polyhedral domains. This problem being piecewisehas some level of natural parallelization and we could further partition thedomain depending on the hardware available. Parallelisation of the codewould require more than just using parallel for statements. We would need toimplement reduction and smartly distribute data in order to get an efficientalgorithm.9.1.5 PreprocessingWe could preprocess the original problem to get different divisions inorder to solve the problem faster. The preprocessing step might vary de-pending on whether we use parallel methods.9.1.6 Overall DomainThe current examples and tests were run when the overall domain wasconvex. We need to run tests when the overall domain is not convex.The current code has been implemented in a way that can be extended toincorporate the functionality that has been mentioned as future possibilities.919.1. FUTURE WORK(a) r1,j(b) s1,1,k (c) s1,3,k(d) skFigure 9.1: Domain of a conjugate illustrating the potential for improve-ment. The conjugates corresponding to vertices and edges that are not apart of the division when computing the convex envelope of each piece ap-pear exactly the same in the final conjugate.92Bibliography[AB10] Kurt Anstreicher and Samuel Burer. “Computable representa-tions for convex hulls of low-dimensional quadratic forms”.Math-ematical Programming 124 (2010), pages 33–43.[AF83] Faiz Al-Khayyal and James Falk. “Jointly constrained biconvexprogramming”. Mathematics of Operations Research 8 (1983).[Ans12] Kurt Anstreicher. “On convex relaxations for quadratically con-strained quadratic programming”. Mathematical Programming136 (2012), pages 233–251.[Bau+08] Heinz Bauschke, Rafal Goebel, Yves Lucet, and Xianfu Wang.“The proximal average: basic theory”. SIAM Journal on Opti-mization (2008).[BC17] Heinz Bauschke and Patrick Combettes. Convex Analysis andMonotone Operator Theory in Hilbert Spaces. Second. 2017.[BLT08] Heinz Bauschke, Yves Lucet, and Michael Trienis. “How to trans-form one convex function continuously into another”. SIAM Re-view (2008).[BM23] Heinz Bauschke and Walaa Moursi. An Introduction to Convex-ity, Optimization, and Algorithms. SIAM, 2023.[Bre89] Yann Brenier. “Un algorithme rapide pour le calcul de trans-forme´es de Legendre-Fenchel discre`tes”. Comptes Rendus de l’Aca-de´mie des Sciences. Se´rie 1, Mathe´matique (1989).[Cor96] Lucilla Corrias. “Fast Legendre-Fenchel Transform and Appli-cations to Hamilton-Jacobi Equations and Conservation Laws”.SIAM Journal on Numerical Analysis 33 (1996).[Cra89] Yves Crama. “Recognition problems for special classes of poly-nomials in 0–1 variables”.Mathematical Programming 44 (1989),pages 139–155.93BIBLIOGRAPHY[GKL14] Bryan Gardiner, Jakee Khan, and Yves Lucet. “Computing thepartial conjugate of convex piecewise linear-quadratic bivari-ate functions”. Computational Optimization and Applications(2014).[GL11] Bryan Gardiner and Yves Lucet. “Graph-Matrix Calculus forComputational Convex Analysis”. In: Fixed-Point Algorithmsfor Inverse Problems in Science and Engineering. Edited byHeinz H. Bauschke, Regina S. Burachik, Patrick L. Combettes,Veit Elser, D. Russell Luke, and Henry Wolkowicz. Springer NewYork, 2011, pages 243–259. url: https://doi.org/10.1007/978-1-4419-9569-8_12.[GL13] Bryan Gardiner and Yves Lucet. “Computing the conjugate ofconvex piecewise linear-quadratic bivariate functions”. Mathe-matical Programming 58 (2013), pages 249–272.[HL07] 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 3 (2007), page 657.[HL18] Tasnuva Haque and Yves Lucet. “A linear-time algorithm tocompute the conjugate of convex piecewise linear-quadratic bi-variate functions”. Computational Optimization and Applications70 (2018).[HL93] Jean-Baptiste Hiriart-Urruty and Claude Lemarechal. ConvexAnalysis and Minimization Algorithms II: Advanced Theory andBundle Methods. Volume 306. Grundlehren der mathematischenWissenschaften. Springer, 1993.[KL20] Deepak Kumar and Yves Lucet. “Towards the Biconjugate ofBivariate Piecewise Quadratic Functions”. In: Optimization ofComplex Systems: Theory, Models, Algorithms and Applications.Springer International Publishing, 2020.[Kum19] Deepak Kumar. “Conjugate of some rational functions and con-vex envelope of quadratic functions over a polytope”. Master’sthesis. University of British Columbia, 2019.[LBT09] Yves Lucet, Heinz Bauschke, and Mike Trienis. “The piecewiselinear-quadratic model for computational convex analysis”. Com-putational Optimization and Applications 43.1 (2009), pages 95–118. url: https://doi.org/10.1007/s10589-007-9124-y.94BIBLIOGRAPHY[Lin05] Jeff Linderoth. “A simplicial branch-and-bound algorithm forsolving quadratically constrained quadratic programs”. Mathe-matical Programming 2 (2005), pages 251–282.[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 (2016).[Luc06] Yves Lucet. “Fast Moreau envelope computation”. NumericalAlgorithms 43 (2006), pages 235–249.[Luc10] Yves Lucet. “What Shape Is Your Conjugate? A Survey of Com-putational Convex Analysis and Its Applications”. SIAM Review52 (2010).[Luc17] Yves Lucet. Computational Convex Analysis (CCA) numericallibrary. https://github.com/ylucet/CCA. 2017.[Luc21] Yves Lucet. Computational Convex Analysis for bivariate piece-wise linear-cubic functions (CCA2) numerical library. https://github.com/ylucet/CCA2. 2021.[Luc96] Yves Lucet. “A fast computational algorithm for the Legendre-Fenchel transform”. Computational Optimization and Applica-tions 6 (1996), pages 27–57.[Luc97] Yves Lucet. “Faster than the Fast Legendre Transform, the Linear-time Legendre transform”. Numerical Algorithms 16 (1997).[MK87] Katta Murty and Santosh Kabadi. “Some NP-complete prob-lems in quadratic and nonlinear programming”. MathematicalProgramming 39 (1987), pages 117–129.[Mor65] Jean-Jacques Moreau. “Proximite´ et dualite´ dans un espace hilber-tien”. Bulletin de la socie´te´ mathe´matique de France (1965).[NN94] Yurii Nesterov and Arkadii Nemirovskii. Interior Point Polyno-mial Methods in Convex Programming. 1994.[Roc70] R. Tyrrell Rockafellar. Convex Analysis. Princeton UniversityPress, 1970.[RW98] R. Tyrrell Rockafellar and Roger J-BWets. Variational Analysis.Springer, 1998.95BIBLIOGRAPHY[SA90] Hanif Sherali and Amine Alameddine. “An explicit characteriza-tion of the convex envelope of a bivariate bilinear function overspecial polytopes”. Annals of Operations Research 25 (1990),pages 197–209.[Sah74] Sartaj Sahni. “Computationally related problems”. SIAM Jour-nal on Computing 3 (1974), pages 262–279.96Appendix ACode demonstrationGiven a PLQ function with two pieces, we demonstrate how to create aPLQ function and get the final output.In the following code, we create a bivariate PLQ function with two pieces.First we create symbolic variables x and y in MATLAB and then create afunction f in terms of x and y. Then we input the vertices of the polyhedralregions for the two pieces and store them in d(1) and d(2). Now that wehave the function and domain, we create two pieces p(1) and p(2) and finallythe plq function and store it in variable - example1.x = sym( ’ x ’ )y = sym( ’ y ’ )f=symbol icFunct ion (x∗y ) ;d(1)=domain ([−5 ,−4;0 ,−4;1 ,3 ;−5 ,5] , x , y ) ;d(2)=domain ( [ 0 , − 4 ; 2 , 0 ; 2 , 1 ; 1 , 3 ] , x , y ) ;p (1 ) = p l q 1p i e c e (d ( 1 ) , f ) ;p (2 ) = p l q 1p i e c e (d ( 2 ) , f ) ;example1 = plq (p ) ;Once we have the input, we can invoke the conjugate function to get theconjugate of the function. The input and output variables are the same asFigure A.1: Function f(x, y) = xy with piecewise polyhedral domain.97Appendix A. Code demonstrationwe are updating fields inside the class.example1 = example1 . conjugateNow we need the output to interpret our results.example1 . printDomainMaple ;example1 . p r i n t ;example1 . pr intLatex ;Line 1 outputs the domain of the convex envelopes, conjugates and max-imum conjugate in a format to generate diagrams in maple. The Figurescreated using this command are displayed in chapters (3-7).Line 2 outputs the convex envelope, conjugate and maximum in theMATLAB command window.Line 3 outputs the convex envelope, conjugate and maximum in LATEXformat.98Appendix BSymbolic verification forproofs in Chapter 7.For Proposition 7.6, we have symbolically verified the expressions inMATLAB. Code is given below. In order to verify “delta” in the codeevaluates to 0.m = sym( ’m’ )q = sym( ’ q ’ )ψ01 = sym(′ψ′01);ψ02 = sym(′ψ′02);ψ03 = sym(′ψ′03);ψ11 = sym(′ψ′11);ψ12 = sym(′ψ′12);ψ13 = sym(′ψ′13);ψ21 = sym(′ψ′21);ψ22 = sym(′ψ′22);ψ23 = sym(′ψ′23);t0 = (−ψ01 −mψ02)/(2(ψ11 +mψ12));t1 = 1/(2(ψ11 +mψ12));t2 = m/(2(ψ11 +mψ12));γ10 = t1 ∗ (ψ23 + qψ22)/(ψ11 +mψ12);γ01 = t2 ∗ (ψ23 + qψ22)/(ψ11 +mψ12);γ00 = (t0 ∗ (ψ23 + qψ22)− ψ13 − q(ψ12))/(ψ11 +mψ12);ζ11 = −(ψ11γ10 +mψ12γ10)2/(ψ23 + qψ22) + γ10;ζ12 = −(2(ψ11γ01+mψ12γ01)(ψ11γ10+mψ12γ10))/(ψ23+ qψ22)+γ01+mγ10ζ22 = −(ψ11 ∗ γ01 +mψ12γ01)2/(ψ23 + qψ22) + γ01mζ10 = −2(ψ11γ01 +mψ12γ10)(ψ13 + ψ11γ00 + ψ12(q +mγ00))/(ψ23 + qψ22)−mψ02γ10 + γ00 − ψ01γ10;99Appendix B. Symbolic verification for proofs in Chapter 7.ζ01 = −(2(ψ11γ01+mψ12γ01)(ψ13+ψ11γ00+ψ12(q+mγ00)))/((ψ23+qψ22))−mψ02γ01 − ψ01γ01 +mγ00 + q;ζ00 = −(ψ13+ψ11γ00+ψ12(q+mγ00))2/(ψ23+qψ22)−ψ03−ψ01γ00−ψ02(q+mγ00);D = ζ212 − 4ζ11ζ22s imp l i f yF ra c t i on (D)s1 = sym(′s′1)s2 = sym(′s′2)d i sp ( ’ quad ’ )fq = simplifyFraction(ζ11s21 + ζ12s1s2 + ζ22s22 + ζ10s1 + ζ01s2 + ζ00)[cx, tx] = coeffs(fq, [s1, s2])δ = simplifyFraction(cx(2)2 − 4cx(1)cx(4))For Proposition 7.8, symbolic code for verification is given below.m = sym( ’m’ )q = sym( ’ q ’ )ψ01 = sym(′ψ′01);ψ02 = sym(′ψ′02);ψ03 = sym(′ψ′03);ψ11 = sym(′ψ′11);ψ12 = sym(′ψ′12);ψ13 = sym(′ψ′13);ψ21 = sym(′ψ′21);ψ22 = sym(′ψ′22);ψ23 = sym(′ψ′23);s1 = sym(′s′1)s2 = sym(′s′2)% Frac t i ona l partζ00 = (ψ21 +mψ22)2 ;δ1 = −2(ψ13ψ21−ψ11ψ23+mψ13ψ22−mψ12ψ23−qψ11ψ22+qψ12ψ21)(ψ01ψ21+ψ211 +m(ψ212m+ ψ01ψ22 + ψ02ψ21 + 2ψ11ψ12 + ψ02ψ22m));100Appendix B. Symbolic verification for proofs in Chapter 7.δ0 = 2ψ11ψ23(ψ11 +mψ12)− 2ψ13ψ21(ψ11 +mψ12)− ψ03ψ21(ψ21 +mψ22) +ψ01ψ23(ψ21 +mψ22) − 2mψ13ψ22(ψ11 +mψ12) + 2mψ12ψ23(ψ11 +mψ12) −mψ03ψ22(ψ21 +mψ22) +mψ02ψ23(ψ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);si1 = 2 ∗ (ψ21 +mψ22)(ψ13ψ21 − ψ11ψ23 +mψ13ψ22 −mψ12ψ23 − qψ11ψ22 +qψ12ψ21)s1 + 2m(mψ22 + ψ21)(ψ13ψ21 − ψ11ψ23 + mψ13ψ22 − mψ12ψ23 −qψ11ψ22 + qψ12ψ21)s2 + δ1;si12 = −(ψ21 + mψ22)s1 − m(ψ21 + mψ22)s2 + ψ01ψ21 + ψ211 + m(mψ212 +ψ01ψ22 + ψ02ψ21 + 2ψ11ψ12 +mψ02ψ22);si0 = (−ψ23(ψ21 +mψ22) − qψ22(ψ21 +mψ22)) ∗ s1 + (qψ21(ψ21 +mψ22) −mψ23(ψ21 +mψ22))s2 + δ0;a = sym( ’ a ’ )b = sym( ’ b ’ )c = sym( ’ c ’ )% Linearfl = as1 + bs2 + c% f r a c t i o n a l form = l i n e a rd i sp ( ’ f rac−l i n ’ )temp = simplifyFraction(((si21/(ζ200si12)))− (fl − si0)2)[cx, tx] = coeffs(temp, [s1, s2])simplifyFraction(cx(2)2 − 4 ∗ cx(1) ∗ cx(4))101Appendix CCode - TestingThe overall code is currently approximately 8,000 lines of MATLAB codewithout tests.C.1 Unit testingWe have unit tested individual methods in the classes along with in-tegration testing where required. For example class region has about fiftymethods, each requiring multiple tests along with integrated testing.C.2 VerificationAs a verification for correctness of piecewise functions, we have checkedthat the functions on adjacent boundaries are the same on the boundaries.Also some derivations like equation of parabolic curve and differentials wereevaluated with different derivations to check that we get the same answers.C.3 Overall testingFor overall testing we have run the same problem with different divisionsas shown in 7.13 and 3.3. We also triangulated the same problem and gotthe same results.102Appendix DDerivations for Step 1D.1 Steps to solve the subproblems.D.1.1 Step (i)We eliminate one variable using the equation ηh(a, b) = ηw(a, b) andreduce the subproblem to a one-dimensional problem on the basis of thefunctional form of the functions ηh(a, b), ηw(a, b) over Shwr .The functions ηh(a, b), ηw(a, b) are linear over Shwr Here ηh(a, b) =f(xvi(a, b), yvi(a, b))−axvi(a, b)−byvi(a, b) and ηw(a, b) = f(xvr(a, b), yvr(a, b))−axvr(a, b)−byvr(a, b) for distinct vertices vi and vr of the polytope, P . Henceusing the equality, we can eliminate variable a and are only left with variableb.The function ηh(a, b) is nonlinear while ηw(a, b) is linear over Shwr Recallsj(a, b) is the unconstrained minimum of fj . For (a, b) ∈ Dh, let zj = a+bmjand definetj(a, b) = fj(sj(a, b))− zjsj(a, b). (D.1)Thustj(a, b) = fj(sj(a, b))− (a+ bmj)sj(a, b).From (5.8) we haveηj(a, b) = f(xj(a, b), yj(a, b))− axj(a, b)− byj(a, b),andy = mjx+ qj .Thereforeηj(a, b) = f(xj(a, b))− axj(a, b)− b(mjxj(a, b) + qj),ηj(a, b) = f(xj(a, b))− (a+ bmj)xj(a, b)− bqj .and we getηj = tj(zj)− bqj .103D.1. Steps to solve the subproblems.Thus the nonlinear equation becomestj(zj)− bqj = f(xvi(a, b), yvi(a, b))− axvi(a, b)− byvi(a, b).Substituting for a, using zj = a+ bmj , we gettj(zj)− bqj = f(xvi , yvi)− (zj − bmj)xvi − byvi .We now have an equation in only one variable b.Functions ηh(a, b) and ηw(a, b) are nonlinear over Shwr and mh = mw. Inthis case Shwr = Dh ∩ Dw and the edges corresponding to indices h and ware parallel. Thus after substituting a with zh = a + bmh and referring to(D.1) we getth(zh)− bqh = tw(zh)− bqw. (D.2)through which we can eliminate b.Functions ηh(a, b) and ηw(a, b) are nonlinear over Shwr and mh ̸= mw.We have Shwr = Dh ∩Dw, zh = a+ bmh, and zw = a+ bmw. We getb =zh − zwmh −mw . (D.3)With respect to these variables the region becomes a box, Ih ∩ Iw whereIj = {zj : f ′j(xvj1) ≤ f ′j(xvj2)}. Next we defineϕj(zj) = tj(zj)− qh − qwmh −mw zj .Using th(zh)− bqh = tw(zw)− bqw and (D.3) we getth(zh)− tw(zw) = b(qh − qw),th(zh)− tw(zw) = zh − zwmh −mw (qh − qw) = zh(qh − qw)mh −mw − zw(qh − qw)mh −mw ,th(zh)− zh (qh − qw)mh −mw = tw(zw)− zw(qh − qw)mh −mw ,ϕh(zh) = ϕw(zw).The above equation allows us to eliminate one variable as it can be provedthat ϕj is invertible over Ij . Thuszw = ϕ−1w (ϕh(zh)). (D.4)104D.1. Steps to solve the subproblems.D.1.2 Step (ii)After eliminating one variable as shown in the earlier step, we can showthat the feasible region of each subproblem can be reduced into a finite unionof intervals with respect to the remaining variable. This means that eachsubproblem given by Equation (5.11), can be split into a finite number ofone-dimensional problems with bound constraints.We discuss the case Shwr = Dh∩Dw and mh ̸= mw as the other cases aresimpler and lead to the same conclusion. In this case we replace the originalvariables a and b with zh and zw and then eliminate zw. With respect tothe new variables, Shwr is defined asf ′h(vh1) ≤ zh ≤ f ′h(vh2),f ′w(vw1) ≤ zw ≤ f ′w(vw2).After substituting for zw using (D.4), we combine the above two inequalitiesto getϕ−1h (ϕw(f′w(vw1))) ≤ zh ≤ ϕ−1h (ϕw(f ′w(vw2))).Thus the set Shwr turns out to be the simple interval with respect to zh. Itis less obvious to show that all constraints ηh(a, b) ≤ ηk(a, b),∀k ̸= h,w cor-respond to a finite union of intervals with respect to the remaining variable.To this aim the system: ηh(a, b) = ηw(a, b),ηh(a, b) ≤ ηk(a, b),(a, b) ∈ Shwr .corresponds to a finite union of intervals with respect to the remaining vari-able. The extreme endpoints of these subintervals correspond to those al-ready shown on the border of Shwr are solutions to{ηh(a, b) = ηw(a, b),ηh(a, b) = ηk(a, b).We now need to prove that the number of solutions to this system is finite.Then we can repeat this procedure for all the inequalities ηh(a, b) ≤ ηk(a, b).For each we get a finite union of subintervals for the single remaining variableand hence the subproblem either has no solution or is reduced to a finitecollection of subintervals with respect to the remaining variable.105D.1. Steps to solve the subproblems.D.1.3 Step (iii)In this section, we discuss the solution of subproblem for different func-tional forms of the functions ηh and ηw.Case ηh, ηw linear In this case we were left with a finite number ofsubproblems of the formmaxdl≤b≤du{ξ1(x, y)b+ ξ0(x, y)}, (D.5)where ξ1 and ξ0 are linear functions with respect to x and y. So we getthe optimal solution as b∗ = dl if ξ1(x, y) ≤ 0 and b∗ = du otherwise.The optimal value is a linear function and we have proved the followingobservation.Observation D.1. When both ηh and ηw are linear, (D.5) has an optimalvalue which is a piecewise linear function with respect to x and y.Case ηh nonlinear, ηw linearIn this case the equationηh(a, b) = ηw(a, b).is equivalent to (D.1). After substitution we are left with the sub-problemmaxdl≤zh≤duξ2(x, y)th(zh) + ξ1(x, y)zh + ξ0(x, y), (D.6)where ξ0, ξ1 and ξ2 are linear functions with respect to x and y. If ξ2 ≤ 0,the optimal solution is either dl or du and the optimal value is a linearfunction. Otherwise, we can get the maximum of the objective by settingthe derivative to zero. In this case we get the optimal value asξ2(x, y)fh(ξ1(x, y)ξ2(x, y))+ ξ0(x, y). (D.7)From this we get the following observation.Observation D.2. When ηh is nonlinear and ηw is linear, there exists apolyhedral subdivision of P such that the optimal value of the sub-problem (D.6)is either a linear function with respect to x, y, or has the functional formgiven in equation (D.7).Case ηh and ηw nonlinear, mh = mw106D.1. Steps to solve the subproblems.In this case the edges are parallel and after eliminating b using (D.2) weget the subproblemmaxdl≤zh≤duth(zh) + zh(x) +th(zh)− tw(zh)qh − qw (y −mhx− qh). (D.8)Now setting λ(x, y) = qh+mhx−yqh−qw ∈ [0, 1],∀(x, y) ∈ P , (D.8) can be rewrittenasmaxdl≤zh≤duλ(x, y)tw(zh) + [(1− λ(x, y)th(zh))] + zhx. (D.9)As this equation is strictly concave, we can solve it by setting the firstderivative to zero. Now if that has no solution or the solution lies outsidethe interval [dl, du], then the optimal solution of the subproblem is either dlor du and the optimal value is a linear function. Otherwise we get the uniquesolutions as z∗h(x, y) where dl ≤ z∗h(x, y) ≤ du. Thus we are interested in(x, y) ∈ P where the first derivative vanishes.In view of1. λ(x, y) ∈ [0, 1] , ∀(x, y) ∈ P , and2. t′h, t′w strictly decreasing.We have this set equal to the polytope{(x, y) ∈ P∣∣∣∣{λ(x, y)t′w(dl) + (1− λ(x, y))t′h(dl) + x ≥ 0λ(x, y)t′w(du) + (1− λ(x, y))t′h(du) + x ≤ 0}.The optimal value, defined on the polytope is given byλ(x, y)t′w(z∗h(x, y)) + (1− λ(x, y))t′h(z∗h(x, y)) + xz∗h(x, y),where z∗h(x, y) is implicitly defined.Case ηh and ηw non-linear, mh ̸= mwWhile solving this case we go to the reformulation in original variablesa and b recalling the definitions for ηh and ηw, we getmaxa,b,c{fh(sh(a, b))− (a+ bmh)sh(a, b)− bqh + ax+ by}subject tofh(sh(a, b))− (a+ bmh)sh(a, b)− bqh + ax+ by =fw(sw(a, b))− (a+ bmw)sw(a, b)− bqw + ax+ bym107D.1. Steps to solve the subproblems.τ1 ≤ sh(a, b) ≤ τ2,τ1 ≤ sh(a, b) ≤ τ2.We then solve the subproblems using KKT conditions and obtain thefollowing conclusion. The intersection of these linear inequalities with Pdefines a polytope over which an interior optimal solution and the corre-sponding optimal value of the subproblem are implicitly defined.D.1.4 Step (iv)The convex envelope is the maximum of all the solutions found in theearlier step. It has been shown that the convex envelope function is charac-terized by a polyhedral division such that it has a specific functional formover each member of the subdivision. In case of the bilinear form we haveexplicit functional forms.108
Thesis/Dissertation
eng
Computer Science
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Lucet, Yves
Computing the conjugate of nonconvex bivariate piecewise linear-quadratic functions
Text