An Approximation Algorithm for Optimal Piecewise Linear Interpolations of Bounded Variable Products

We investigate the optimal piecewise linear interpolation of the bivariate product xy over rectangular domains. More precisely, our aim is to minimize the number of simplices in the triangulation underlying the interpolation, while respecting a prescribed approximation error. First, we show how to construct optimal triangulations consisting of up to five simplices. Using these as building blocks, we construct a triangulation scheme called crossing swords that requires at most - times the number of simplices in any optimal triangulation. In other words, we derive an approximation algorithm for the optimal triangulation problem. We also show that crossing swords yields optimal triangulations in the case that each simplex has at least one axis-parallel edge. Furthermore, we present approximation guarantees for other well-known triangulation schemes, namely for the red refinement and longest-edge bisection strategies as well as for a generalized version of K1-triangulations. Thereby, we are able to show that our novel approach dominates previous triangulation schemes from the literature, which is underlined by illustrative numerical examples.


Introduction
We consider optimal piecewise linear (pwl.)interpolations of the bilinear nonconvex function F : R 2 → R, F(x, y) = x y over the rectangular domain D = [ x, x] × [ ȳ, ȳ].A pwl. interpolation f : D → R of F is uniquely defined by its underlying triangulation T := {T 1 , . . ., T k } ⊆ R 2 , k ∈ N of the domain D, where f is linear over each simplex T i .The approximation error between F and f is defined as the maximum pointwise deviation over D. For any prescribed approximation accuracy ε > 0, we are interested in triangulations that lead to pwl.interpolations with an approximation error less than ε and that are minimal with respect to the number of simplices they contain.To the best of our knowledge, finding triangulations which are optimal in this sense is still an open problem.Furthermore, there have been no attempts so far to prove approximation guarantees for existing triangulation schemes in the literature either.

Contribution.
In this paper, we make several important steps toward efficient triangulation schemes for the pwl.interpolation of bilinear functions.At first, we construct triangulations that provably minimize the approximation error with up to five simplices.We show that these triangulations are in turn optimal for approximation accuracies up to ( √ 5−2) /4 ≈ 0.059 over the unit box.Using these basic triangulations as building blocks, we derive an approximation algorithm for the optimal triangulation problem which we call crossing swords, based on its underlying geometric idea.We prove that this scheme can be used to create triangulations consisting of at most √ 5 /2times the number of simplices in an optimal triangulation.To be more precise, the crossing swords algorithm is an √ 5 /2approximation algorithm for all ε i = 1 /16i with i ∈ N.For any intermediate value, the approximation guarantee is √ 5 /2 + 4 √ 5ε.In this sense, the crossing swords scheme is an asymptotic √ 5 /2-approximation algorithm for small prescribed accuracies.Furthermore, we prove that crossing swords produces optimal triangulations for all ε i = 1 /16i with i ∈ N if we require each simplex to have one axis-parallel edge.We will show that crossing swords triangulations are superior to the most-widely used triangulation schemes in the literature, i.e., K1, J1 [21], longestedge bisection [10,14], maximum error bisection [13,20] and red refinement [4,8].Although we show that longest-edge bisection is also a √ 5 /2approximation algorithm, it produces more simplices than crossing swords expect for certain specific values of ε, where the two are equivalent.Furthermore, we introduce a generalized version of the K1-triangulation scheme and show that it is a √ 5-approximation algorithm.This also applies to the red refinement scheme.Since the approximation accuracies in crossing swords triangulations can be adjusted more finely, our scheme is by a factor of two better than previous methods in many cases.The overall dominance of the crossing swords scheme is underlined by numerical results for an indicative example.
Literature Review.Pwl.interpolations of bilinear functions are studied in the literature in various contexts.On the one hand, such approximations are important in many applications.For example, in computer visualization, nonlinear surfaces are usually represented by pwl.functions based on triangulations.This also applies to nonlinear shapes in architectural design, which in some cases can be realized in practice only via pwl.approximations, cf.[19].On the other hand, several very efficient algorithmic frameworks for solving quadratic optimization problems rely on pwl.approximations [1,6,8,13,14].However, finding optimal pwl.interpolations of x y over rectangular domains is still an unsolved problem.The only attempt to find provably minimal triangulations is represented by [15].Here, the author develops a mixed-integer quadratic program which models the optimal triangulation problem.However, due to its size and inherent complexity, the presented model cannot even be solved for trivial instances by state-of-the-art solvers.Besides that, there exist several heuristic triangulation approaches in the literature.In [23], the author uses regular uniform triangulations, namely J1-and K1-triangulations, which were first presented in [21].In [11] and [20], the authors go a step further and develop triangulations which are specifically designed for the approximation of x y over rectangular domains.In [9,13,14], triangulations are constructed adaptively, using red, maximum error or longest-edge refinements.For all mentioned triangulation schemes, we derive relations between the approximation accuracy and the resulting number of simplices.Most importantly, we introduce a new, more efficient triangulation scheme for which we can prove that it dominates the previously known methods.
On a more general basis, the authors of [19] are interested in optimal triangulations of x y over the whole plane R 2 .They show how to construct so-called optimal simplices, i.e., triangles which fulfill a prescribed approximation accuracy while maximizing their area.They prove that R 2 can be tiled with optimal simplices only and thus obtain an optimal triangulation.A natural idea would be to use optimal simplices to triangulate rectangular domains as well.However, due to their geometry, it is not possible to use solely optimal simplices to triangulate rectangular regions.Instead, we will use this idea to derive a lower bound on minimal triangulations over any given compact domain.This bound is essential in proving the above-mentioned approximation guarantees for the different triangulations schemes.

Structure.
Our work is structured as follows.In Sect.2, we introduce the general notation and concepts that are used throughout the article.After that, in Sect. 3 we focus on optimal triangulations for the pwl.interpolation of x y over rectangular domains.We show that the problem can be reduced from general rectangular domains to the unit box.In Sect.3.1, we construct optimal triangulations consisting of up to five simplices.We use them as building blocks to develop our √ 5 /2approximation algorithm called crossing swords in Sect.3.2.It even yields optimal triangulations if we require that each simplex has one axis-parallel edge.In addition, we also prove approximation guarantees for a generalized version of K1-triangulations, the red refinement and the longest-edge bisection scheme.In Sect.3.4, we prove the dominance of the crossing swords triangulation compared to the known triangulations from the literature and provide numerical results on an exemplary instance that underline the theoretical results.Finally, we conclude in Sect. 4.

Piecewise Linear Functions and Approximations
We start by introducing the basic concepts of triangulations and piecewise linear functions and discuss their use in approximating nonlinear functions.For the sake of simplicity, we restrict ourselves to continuous functions over polytopal domains.
A function is called piecewise linear (pwl.) if it is linear over each element of a given domain partition.To partition a domain, it is possible to use any family of polytopes.However, in practice most often triangulations are used, see, e.g., [22].This is without loss of generality, since a pwl.function defined with respect to a partition by polytopes can always be represented by a pwl.function over a triangulation, namely by triangulating each polytope.Therefore, we define pwl.functions over triangulations.In the following, we formally introduce the relevant definitions in this context.Throughout this work, we use the notation V (P) for the vertex set and A(P) for the area of a polytope P ⊂ R d .The following definitions are formulated in a general way for pwl.function in R d .In Sect.3, we only consider the special case d = 2 and therein triangulations formed by full-dimensional simplices, i.e., the special case A triangulation is a partition consisting of full-dimensional simplices.

Definition 2.2 A finite set of full-dimensional simplices
i=1 T i and the intersection T i ∩ T j of any two simplices T i , T j ∈ T is a proper face of T i and T j .Further, we denote the nodes N (T ) of the triangulation T as N (T ) We use the concept of triangulations to define pwl.functions.
Pwl. functions can be used to approximate nonlinear functions.We consider the case where the pwl.approximation coincides with the nonlinear function at the nodes of the triangulations.In this way, we can ensure continuity of the approximation if the nonlinear function is continuous.

Definition 2.4
Let P ⊂ R d be a polytope, and let T := {T 1 , . . ., T k }, k ∈ N, be a triangulation of P. We call a pwl.function g : P → R a pwl.interpolation of a continuous function Usually, the error of a pwl.approximation is measured by the maximum absolute pointwise deviation between the pwl.approximation itself and the nonlinear function to be approximated; see, e.g., [11,13,18,23].In the following, we also use this definition of the approximation error.
Given some ε > 0, we call g an ε-interpolation and T an ε-triangulation if the approximation error is smaller than or equal to ε.
Finally, we formulate the concept of optimal triangulations for pwl.interpolations of nonlinear functions; see also [15].Definition 2.6 Let P ⊆ R d be a polytope and g : P → R a pwl.ε-interpolation of G : P → R with respect to the triangulation T .We call T ε-optimal if |T | is minimal among all ε-triangulations.
It is not obvious how to determine ε-optimal triangulations in general.We tackle this problem for bilinear functions in the following.

Triangulations for the Interpolation of xy Over Box Domains
In this section, we focus on finding optimal triangulations for pwl.interpolations of the bilinear function i.e., we treat the following problem: Problem 1 Given some ε > 0, find an ε-optimal triangulation of D w.r.t.F.
Our main contribution will be the derivation of a novel approximation algorithm for Problem 1.We begin by stating a known result from the literature for the approximation of F over a single simplex.Lemma 3.1 (Approximation error) [19] Given a pwl.interpolation f : D → R of F defined by a triangulation T of D, the approximation error ε f ,F (T ) over a simplex T ∈ T is attained at the center of one of its facets.Further, if (x 0 , y 0 ) and (x 1 , y 1 ) are the endpoints of a facet, we have Fig. 1 Geometrical visualization of the approximation error on an edge of a simplex Since the error on each simplex is attained on one of its facets, we know the following for the approximation error over D.
Lemma 3.2 Let f : D → R be a pwl.interpolation of F defined by a triangulation T of D. Then, the approximation error ε f ,F (T ) is attained on a facet of some simplex T of T that is not on the boundary of D.
Proof The proof follows directly from Lemma 3.1 and the fact that all boundary facets are axis-parallel.This implies that the approximation error is zero there as either x 1 = x 0 or y 1 = y 0 holds in Eq. ( 1).
A geometric view on Lemma 3.1 is sketched in Fig. 1.Here, the maximum approximation error on the facet e v 1 ,v 2 is illustrated.The yellow rectangle is the facet enclosing axis-parallel rectangle, and the area of the red rectangle equals the maximum approximation error on e v 1 ,v 2 .It is one quarter of the area of the yellow rectangle.We use Lemma

be the triangulation of the box D such that each simplex T D ∈ T D corresponds to a simplex T U ∈ T U via the mapping L on the respective vertex sets, i.e., V (T
. The endpoints of e D are nodes of T D and therefore the result of a linear mapping of two vertices (x U 0 , y U 0 ) and (x U 1 , y U 1 ) that are the endpoints of some face e U in T U : Further, let ε e ∈ [0, ε] be the approximation error over e U .The approximation error over e D , attained at its center (x * , y * ), is calculated as follows: Thus, the approximation error on each facet of T D is ν times the error on the corresponding facet in T U .Since we assumed that T U is an ε-interpolation, ε e ν ≤ εν holds for each facet e of T D .This means T D is an νε-triangulation.
Next, we prove that if T U is ε-optimal, then T D is νε-optimal.Assume T U is εoptimal and there exists a νε-triangulation If we apply the inverse of L to the node set of S D , we obtain a triangulation S U of U , i.e., for every simplex . The endpoints of e U are nodes of S U and therefore the result of a linear mapping of two vertices that are the endpoints of some face e D in S D : Further, let νε e with ε e ∈ [0, νε] be the approximation error over e D .The approximation error over e U , attained at its center (x * , y * ), is calculated as follows: Thus, the approximation error on each facet of S U is less than or equal to ε.This means that S U is an ε-triangulation of U with |S U | < |T U |, which contradicts the assumption that T U is ε-optimal.It follows that under the assumption that T U is an ε-optimal triangulation, there is no νε-triangulation of D with fewer simplices than T D , proving that T D is νε-optimal.
Note that in Lemma 3.3 the bound νε is tight if ε is a tight bound for the approximation over U , i.e., if the approximation error ε is attained at some point (x, y) ∈ U , then an approximation error of νε is attained at L(x, y) ∈ D. As a result of Lemma 3.3, we can consider the approximation of F over U w.l.o.g. in the following: Problem 2 Given some ε > 0, find an ε-optimal triangulation of U w.r.t.F.

Solving the Optimal Triangulation Problem for up to Five Simplices and a General Lower Bound
In the following, we solve Problem 2 for a fixed number of simplices up to five and give a general lower bound with respect to the approximation quality ε > 0. Finding a general scheme that solves Problem 2 for arbitrary values of ε is an open problem.However, in [15] it is shown that Problem 2 can be formulated as a mixed-integer quadratically constrained program (MIQCP).The author takes advantage of the fact that the approximation error is attained on one of the facets of the simplices and can therefore exploit the representability of a triangulation by a fully connected planar graph.This representation allows for the modeling of Problem 2 by an MIQCP with a finite number of constraints and variables.To the best of our knowledge, this is the only work trying to determine provably ε-optimal triangulations over box domains.Unfortunately, due to the size of the resulting MIQCP, this approach is computationally intractable even for trivial instances.Nevertheless, we will later use its underlying idea of representing a triangulation as a planar graph to prove the optimality of triangulations consisting of up to five simplices.
A lower bound for ε-optimal triangulations.We begin our examination of Problem 2 in deriving a general lower bound for ε-optimal triangulations.This lower bound enables us later to determine approximation guarantees for specific triangulation schemes.The idea of the presented lower bound goes back to the work of Pottmann et al. [19], who studied optimal triangulations of the plane R2 .The authors showed that an ε-optimal triangulation of R 2 , in the sense that the individual simplices have maximal area, is obtained by using so-called ε-optimal simplices.An ε-optimal simplex satisfies a prescribed accuracy ε while maximizing its area.The area of an ε-optimal simplex is 2 √ 5ε.As it is possible to tile the R 2 with ε-optimal simplices only, such triangulations are optimal.In Fig. 2, we illustrate two different 0.25-optimal simplices.For more information on ε-optimal simplices, we refer the reader to [19] and [2].We use the following lemma from [5], to give a general lower bound for the triangulation of polytopal domains.Lemma 3.4 (Lower bound) [5] Let P ⊂ R 2 be a polytopal domain with an area of A(P).Then, any εtriangulation T P of P w.r.t.F contains at least A(P) The proof is straightforward.If we assume that we can triangulate P solely by εoptimal simplices, which have an area of 2 √ 5ε, we obtain the indicated lower bound immediately.For the approximation of F over U , Lemma 3.4 thus sets the following lower bound: Corollary 3.5 An ε-optimal triangulation T of U w.r.t.F requires at least Fig. 2 Optimal simplices with an approximation error of 0.25 However, it is unclear whether or how to use ε-optimal simplices to triangulate U .Furthermore, the lower bound from Lemma 3.4 is not always tight.In [17], the author proves that one cannot triangulate a rectangle with an odd number of simplices such that all simplices have the same area.This means for those prescribed approximation accuracies ε > 0 for which the lower bound is an odd number, we need at least one additional simplex.We now show that this lower bound can even be improved if we require that each simplex has at least one axis-parallel edge.This lemma has also been presented in the dissertation of the fourth author [15].

Lemma 3.6 Let ε > 0 be a prescribed approximation accuracy for the interpolation of F over a simplex T . If T has at least one axis-parallel edge, then the maximum area of T is 4ε.
Proof Let T be a simplex with vertices v 1 = (x 1 , y 1 ), v 2 = (x 2 , y 2 ) and v 3 = (x 3 , y 3 ).W.l.o.g., we assume e v 1 ,v 2 is parallel to the x-axis.The area of T can then be calculated using the areas of the axis-parallel edge-enclosing rectangles of e v 1 ,v 3 and e v 2 ,v 3 : Further, for the edges e v 1 ,v 3 and e v 2 ,v 3 , the approximation error has to be less than or equal to ε, i.e., Now, obviously A(T ) attains its maximum possible value of 4ε if holds.Since we assumed e v 1 ,v 2 to be parallel to the x-axis, y 1 = y 2 holds and we get a simplex of maximum area 4ε if the vertices have the following positions: Corollary 3.7 An ε-optimal triangulation T of U w.r.t.F where each simplex has at least one axis-parallel edge requires 1 4ε simplices.
Euler relations for triangulations.The following optimality proofs exploit the representation of a triangulation by a fully connected planar graph.In this sense, the vertices and facets of the simplices become the nodes and edges of a graph.This allows us to translate the well-known Euler relations for planar graphs to triangulations, see also [3].
Let T be a triangulation of some polytope P ⊂ R 2 , then it holds that and where K is the number of vertices that lie on the boundary of P and E is the number of edges of T .For K , we further know that which together with (2) gives If P is a rectangle, then it obviously also holds that Optimal triangulations with up to five simplices.We use the Euler relations to derive optimal triangulations of U for the interpolation of F under the restriction that the number of simplices is fixed a priori.These triangulations in turn yield ε-optimal triangulations in the sense of Problem 2 for specific value ranges of ε.Let T * n be a triangulation with a minimal approximation error among all triangulations of U consisting of n simplices, and let ε * n be the corresponding minimal approximation error.This means that T * n is an ε-optimal triangulation for all ε ∈ (ε * n+1 , ε * n ].In the following, we determine optimal triangulations T of U for |T | = 2, 3, 4, 5.In both cases, the following linear program defines the optimal position of the fifth node:

We start with
In Problem 7 and the following optimization problems, we do not give specific names to variables.Instead, we use the variable notation x i as it is conform to the standard definition of quadratic programs, which we use in Lemma 3.8.Nevertheless, in Problem 7 and in the subsequent optimization problems, 1  4 x 1 always models the bound on the approximation error that is minimized, see Lemma 3.1 for the calculation of the approximation error.All further variables x i model the axis-lengths of the inner edges which can be chosen freely in the respective configuration.Each constraint models the approximation error on an inner edge of the considered triangulation.The approximation error on an outer edge is always zero.Figures 3(a)-3(i) illustrate which constraint models which edge.In Fig. 3b, we have two inner edges that both have a y-axis length of 1 and an x-axis length of x 2 and 1 − x 2 , resulting in maximum approximation errors on these edges of 1  4 x 2 and 1 4 (1 − x 2 ).Since both errors must be less than or equal to the bound 1  4 x 1 , the constraints in Problem 7 are obtained after transformation: The optimal solution to Problem 7 is x * = 1 2 , 1 2 with a corresponding approximation error of ε * 3 = 1 8 .We refer to these equivalent triangulation by T * 3 , and an example is shown in Fig. 3b.
We continue with |T | = 4: Again, from relations ( 2)-( 6) we know that |N (T )| ∈ {4, 5, 6}.If |N (T )| = 4, we obtain K = 2, which is a contradiction to (6).If |N (T )| = 5, we obtain K = 4 and E = 8.Triangulations that fulfill this property have one inner vertex that is connected to all four corners of U .We can find the optimal position of the inner vertex by solving this quadratic program: In Fig. 3c, we have four inner edges with maximum approximation errors of 1  4 x 2 x 3 , ).All errors have to be less than or equal to 1  4 x 1 , which after transformation leads to the constraints in Problem 8.The optimal solution to Problem 8 is x * = 1 4 , 1 2 , 1 2 , with an approximation error of ε * 4 = 1 16 .We refer to the corresponding optimal triangulation by T * 4 , see Fig. 3c.Finally, if |N (T )| = 6, we obtain K = 6, which means that all nodes have to be on the boundary of U .Due to symmetry, the remaining two free nodes can either be positioned on adjacent, opposite or both on the same boundary facet.If both nodes are placed on the same boundary facet of U , we have a similar case to |T | = 3 and the minimal possible approximation error is 1  8 .An example of this triangulation is shown in Fig. 3d.If the two nodes are placed on adjacent boundary facets, we can minimize the approximation error by solving the following quadratic program: In Fig. 3d, we have three inner edges with maximum approximation errors of 1  4 x 2 x 3 , 1 4 (1 − x 2 ), 1  4 x 2 (1 − x 3 ), and 1 4 (1 − x 3 ).All errors have to be less than or equal to 1  4 x 1 , which after transformation leads to the constraints in Problem 9.The optimal value of Problem 9 is 1 /8.Thus, this configuration can also not be optimal, see Fig. 3e, for example.If the two additional nodes are placed on opposite boundary facets, we can minimize the approximation error by solving the following linear program: In Fig. 3e, we have three inner edges with maximum approximation errors of 1  4 x 2 , 1 4 (1 − x 3 ), 1 4 (x 2 − x 3 ), and 1 4 (1 − x 3 ).All errors have to be less than or equal to 1 4 x 1 , which after transformation leads to the constraints in Problem 9.The optimal solution to Problem 10 is x = ( 1 3 , 1 3 , 2 3 ), with an approximation error of 1  12 .Thus, this configuration can also not yield an optimal triangulation, see Fig. 3f for an example.Consequently, we know that ε * 4 = 1 /16 holds and that T * 4 is optimal.Finally, we consider the case |T | = 5: From relations (2)-( 6), it follows that |N (T )| ∈ {6, 7}.If |N (T )| = 6, we obtain K = 5.This means we have one free node on the boundary of U and one free inner node.In this case, there are three different ways to connect the nodes such that we have a triangulation with exactly five simplices.It is obvious that two of them are not optimal, see Fig. 3g and h.The approximation error of the third configuration is minimized by the following nonconvex quadratically constrained quadratic program (QCQP), visualized in Fig. 3i: In Fig. 3i, we have five inner edges with maximum approximation errors of 1  4 (1 − x 2 ), ).All errors have to be less than or equal to 1  4 x 1 , which after transformation leads to the constraints in Problem 11.By geometric reasoning, we guess as a candidate for an optimal solution, which entails an approximation error of ε * 5 = . In the following, we prove the optimality of x * .In general, we can prove that a vector is globally optimal for a nonconvex QCQP by checking two conditions which together are sufficient.Lemma 3.8 (Global optimality in nonconvex QCQPs) [16] Consider a quadratic program of the form where all Q i are n × n real symmetric matrices and c i ∈ R n .A vector x * ∈ R n is globally optimal for this problem if there exists a vector and hold.
We use Lemma 3.8 to prove that x * is optimal for Problem 11.

Proposition 3.9 The vector x
) T is globally optimal for the nonconvex quadratic program (11).
Proof First, we show that x * is a locally optimal solution, applying the KKT conditions from (12).This means we have to find a vector λ ∈ R 5 that solves the following system: Note that we can neglect the dual variables for the constraints x ∈ [0, 1] 4 , since they must all be zero for x * anyway due to 0 < x * i < 1.We can easily check that λ * = ( ( , 0, 0, 0) T is feasible for (14), and therefore, (x * , λ * ) is feasible for (12).According to Lemma 3.8, (x * , λ * ) is globally optimal for Problem (11) if is strictly upper triangular, its only eigenvalue is 0. It is therefore a positive semidefinite matrix, and consequently, x * is a globally optimal solution to Problem 11.
As a result of Proposition 3.9, the minimal approximation error of the configuration modeled in Problem 11 is ε * 5 = √ 5−2 4 .This means that the corresponding triangulation T * 5 is optimal for the case of five simplices.

Crossing Swords: A √ 5/2-Approximation Algorithm for Optimal Triangulations
We use the optimal triangulations consisting of up to five simplices from Sect.3.1 to derive our novel triangulation scheme crossing swords stated in Algorithm 1.The algorithm takes as input an integer N and returns a triangulation of U with N simplices that has an approximation error of less than 1 4(N −1) .We will use this approximation accuracy depending on the number of simplices to prove that the crossing swords scheme is a

Proof
In the following, we prove the correctness of the approximation error stated in Eq. ( 15) with respect to the modulus of N by 4. For the sake of simplicity, we speak of T * i -triangles ( T * i ) whenever a rectangle is triangulated in fashion of T * i with i ∈ {2, 3, 4, 5}.Via the scaling formulas established in Lemma 3.3, the approximation error of these subtriangulations on arbitrary rectangles D with an area of A(D) is given as follows (cf.Fig. 3): • A(D).
We now distinguish the four cases of Eq. (15), and examples of these with N = 6, 7, 8, 9 are shown in Fig. 5.To this end, we check that the areas of the rectangles sum up to one and that the approximation error postulated in Eq. ( 15) is fulfilled by the returned triangulation T N cs of Algorithm 1. Case N ≡ 0 mod 4: Each U i has an area of A(U i ) = 4 N , and therefore, holds.The approximation error on each T * 4 -triangle, and thus also the one arising from T N cs , is given as Case N ≡ 1 mod 4: The rectangles U 1 , . . ., U n−1 have an area of for the number of simplices in an ε-optimal triangulation from Proposition 3.5, we can prove the proclaimed approximation guarantee of crossing swords triangulations: This upper estimate of the approximation guarantee of crossing swords triangulations can be strengthened when only considering certain discrete values of approximation accuracies.
Theorem 3.13 Let ε := 1 16i for some i ∈ N.Then, calling Algorithm 1 with input 2 -approximation algorithm for Problem 2. Proof By the choice of ε, we have N ε ∈ N. Further, we know from Eq. ( 15) that N ε ≡ 0 mod 4 holds and that Algorithm 1 called with N ε returns an ε-triangulation T N ε cs .Again, by using the lower bound from Proposition 3.5, we can prove the proclaimed approximation quality of crossing swords triangulations: Next, we show that crossing swords produces optimal triangulations for ε := 1 16i with i ∈ N if we additionally require that each simplex has at least one axis-parallel edge.Theorem 3.14 Let ε := 1 16i for some i ∈ N.Then, calling Algorithm 1 with input N ε := 1 4ε yields an ε-optimal triangulation for Problem 2 if we require that each simplex has at least one axis-parallel edge.

123
Proof In Corollary 3.7, we showed that under the condition that each simplex has an axis-parallel edge, the lower bound for the number of simplices in an ε-optimal triangulation is 1 4ε .Further, we know from Eq. ( 15) that N ε ≡ 0 mod 4 holds and that Algorithm 1 called with N ε returns an ε-triangulation T N ε cs , where Note that in Algorithm 1, we do not explicitly state how to construct the rectangular partition of the domain.In fact, a partition with given areas of the rectangles always exists as we explain now.One way to create such a partition is to line up rectangles along the x-axis that have a full height of one.In this arrangement, it is easy to choose the area of these rectangles as required based on their width, whereas it is not obvious how to do for other arrangements.A corresponding version of Algorithm 1 is stated in the appendix as Algorithm 2. However, this arrangement inherently leads to very skinny simplices for large values of N , see Fig. 4 for an example with N = 16.Skinny simplices in turn can lead to numerical difficulties in evaluating the function values of the pwl.interpolation.Although irrelevant from a theoretical point of view, it would be desirable in practice to partition U into rectangles that are "as square as possible" to avoid skinny simplices.In [7], two problems are studied that address this very question of partitioning "as square as possible".One is called PERI-SUM, which aims at minimizing the sum of all perimeters for a given set of rectangular areas.The second one is PERI-MAX, which has the goal to minimize the largest perimeter over all rectangles for a given set of rectangular areas.Both problems are NP-complete.For certain inputs, namely N = 4 i with i ∈ N, all rectangles have the same area and a partition into squares is possible.The longest-edge bisection scheme described in Sect.3.3 yields exactly such crossing swords triangulations with square partition elements for N = 4 i .For general areas as inputs to PERI-SUM, [12] gives a modeling as an MIQCP.In order to find high-quality solutions fast, the authors also present a polynomial time 3 √ 2 -approximation algorithm.Finally, we make the conjecture that Algorithm 1 indeed yields optimal triangulations for a special sequence of errors converging to zero which we previously used in Theorem 3.14.The following conjecture is equivalent to Theorem 3.14 without the requirement of axis-parallel edges and yields triangulations with N ε ≡ 0 mod 4 many triangles.

Conjecture 3.15 Let ε := 1
16i for some i ∈ N.Then, calling Algorithm 1 with input N ε := 1 4ε yields an ε-optimal triangulation for Problem 2. Remark 3. 16 According to Lemma 3.3, we can apply Algorithm 1 and the lower bound from Lemma 3.4 to general box domains Therefore, the approximation guarantees from Theorems 3.12 and 3.13 are also valid for arbitrary box domains.

Approximation Qualities of Widely Used Triangulation Schemes
In this section, we prove approximation guarantees for several popular triangulation schemes from the literature.For this purpose, we proceed analogously to Theorem 3.13 and consider only discrete values of approximation accuracies converging to zero.Generalized K1-triangulations We start by analyzing the uniform triangulation schemes J1 and K1.Furthermore, we develop a generalized version of the K1, for which we prove an approximation guarantee of √ 5.In Fig. 6a and b, we show a K1and a J1-triangulation consisting of 32 simplices each.Both triangulation schemes are defined over a grid that partitions the domain into axis-parallel squares.Each square is triangulated in fashion of T * 2 .Note that J1 and K1 only differ in the orientation of the diagonals.However, the orientation of the diagonals does not matter for the resulting approximation error.Therefore, both triangulation schemes are equally efficient in terms of the number of simplices.Proposition 3.17 For any given uniform rectangular grid, K1-and J1-triangulations have the same approximation error.
In the way, K1-and J1-triangulations have been used in the literature so far; they require the same uniform partitioning along both axes.Therefore, they offer only limited flexibility to generate good triangulations for most prescribed approximation accuracies.To add somewhat more flexibility, we generalize the K1-triangulation scheme by allowing subdivisions of different granularities for each axis-see also Fig. 6c, where we show a generalized K1-triangulation over a 4 × 2-grid.The approximation error of the generalized K1-triangulation scheme is uniquely determined by the underlying grid.Proposition 3.18 Let L be a uniform grid over U with i intervals on the x-axis and j intervals on the y-axis.Further, let T i j K1 be a generalized K1-triangulation of U with respect to L. Then, 123 Proof Each rectangle of the underlying grid has an area of 1 i j .From Lemma 3.3, we know that the approximation error of a T * 2 triangulations is 1  4 times the area of the triangulated rectangle, which finishes the proof.
We can use the lower bound from Lemma 3.4 to prove an approximation guarantee for generalized K1-triangulations.Theorem 3.19 Let ε := 1 4i j for some i, j ∈ N.Then, the generalized K1-triangulation scheme is an √ 5-approximation algorithm for Problem 2.
Proof Let L be a uniform grid over U with i intervals in x-direction and j intervals in y-direction, and let T i j K1 be the corresponding K1-triangulation.From Proposition 3.18, it follows that ε(T As a result, we can prove the proclaimed approximation quality by using the lower bound from Proposition 3.5: Together with Theorem 3.13, this means that generalized K1-triangulations need twice the number of simplices than crossing swords triangulations to fulfill the same approximation accuracy.This observation is visualized in Fig. 6d, where we show how to convert a J1-triangulation with 32 simplices into a crossing swords triangulation of 16 simplices that has the same approximation accuracy by removing certain edges.

Approximation Quality of Iterative Refinement Schemes Starting with
In the following, we show approximation qualities for triangulations generated by refinement schemes.A refinement scheme specifies how the simplices of a given triangulation are subdivided into smaller simplices.In the context of pwl.approximations, a refinement is performed to reduce the approximation error.In the following, we always assume that the initial triangulation of U is T * 2 .Subsequently, we perform iterative refinements on all simplices where the approximation error is greater than an a prescribed accuracy ε > 0 and repeat the refinements until we end up with an ε-triangulation.Longest-edge bisection.The longest-edge bisection subdivides a simplex by adding a new vertex at the center of its longest edge and connects it with the opposite vertex, see Fig. 7. Algorithm 3 in the appendix gives a detailed version of this procedure.We define T i L as the triangulation resulting from the i-th longest-edge refinement of T 0 L := T * 2 .This strategy also results in a T T T Proof The error on an edge e with endpoints (x 0 , y 0 ) and (x 1 , y 1 ) is always attained at its midpoint and has a value of 1  4 |(x 1 − x 0 )(y 1 − y 0 )|.If e is now bisected, the approximation error on the resulting edges reduces by a factor of 1  4 .As we start with T 0 L = T * 2 , the approximation error in both initial simplices is attained at the midpoint of the diagonal edge and is 1  4 .The longest-edge bisection is now performed on both simplices at the midpoint of the diagonal, see Fig. 7b.The resulting four simplices each have one of the outer edges of U as their longest edge, with a constant approximation error of 0 along that edge.However, the maximum approximation error on each simplex is on one of the new diagonal edges and has a value of 1  16 .If we now perform a longest-edge bisection again, we bisect exclusively on axisparallel edges.This is because adding a new node on an axis-parallel edge, along which the pwl.interpolation corresponds to F anyway, leaves the pwl.interpolation unchanged.The resulting eight simplices shown in Fig. 7c therefore have the same approximation error as the four from the previous iteration.In this way, we continue the refinement and obtain an increase in the number of simplices by a factor of two in each refinement iteration, while the approximation error is reduced by a factor of 1  4 in the odd refinement iterations only.This leads to the fact that the triangulation T i L from the i-th longest-edge refinement consists of 2 i+1 simplices and has an approximation error of 1  2 i+3 .Therefore, we obtain the proclaimed approximation guarantee of √ 5/2 only for the resulting triangulation after odd iterations of the longest-edge refinement strategy: 2 .

Red refinement
The red refinement procedure subdivides a simplex into four smaller simplices by adding a new vertex at the center of each edge and new edges that connect these new vertices, see Fig. 8. Algorithm 4 in the appendix gives a detailed version of Proof By Proposition 3.18, the discrete approximation errors attained exactly by the generalized K1-triangulations are all ε i j := 1 4i j with i, ∈ For these approximation errors, we know that generalized K1 is an √ 5-approximation algorithm and that crossing swords is an √ 5( 1 2 + 4ε i j )approximation algorithm.For i, j = 1, i.e., ε 1,1 = 1 4 , both schemes produce equivalent triangulations consisting of two simplices.For i = 1 and j = 2 or i = 2 and j = 1, i.e., ε 1,2 = ε 2,1 = 1 8 , the generalized K1-triangulation produces four simplices, while crossing swords returns only three simplices.For all other values of i and j, ε i, j < 1 8 holds, and therefore, crossing swords has a guaranteed approximation quality of less than √ 5.

Corollary 3.23 (Crossing swords dominates longest-edge bisection)
For any ε > 0, the crossing swords triangulation scheme produces triangulations that have at most the number of simplices as a triangulation constructed by longest-edge bisection.
Proof By Theorem 3.20, longest-edge bisection has an approximation quality of for all ε i := 1 2 i+3 with i ∈ N.These values are a subset of the approximation errors ε j := 1  16 j , for which the crossing swords triangulation also has an approximation quality of  With a similar proof as for Corollary 3.22, we also obtain that crossing swords outperforms red-refinement.
Corollary 3.24 (Crossing swords dominates red-refinement) For any ε > 0, the crossing swords triangulation scheme produces triangulations that have at most the number of simplices as a triangulation constructed by longest-edge bisection.For ε < 1 4 , the number is strictly lower.
In addition to these theoretical proofs, we also give a numerical example that underlines the presented properties.We consider the domain D = [0, 2] × [2,6] and approximation accuracies of {1, 0.5, 0.25, 0.1, 0.05}.The results are summarized in Table 1.Here, we list the prescribed approximation accuracy ε, the number of simplices and the exact approximation error of the respective triangulation.The stated results underline the efficiency and adaptability of the crossing swords triangulation scheme.For the finest approximation accuracy ε = 0.05, the crossing swords scheme produces a triangulation with only 60 simplices, while the other schemes need a factor of at least 1.5 times more simplices.The next best triangulation scheme is that by [20], which is based on a maximum error bisection refinement strategy.However, they consider general pwl.approximations instead of interpolations.Thus, since they solve a generalization of Problem 2, the results are not directly comparable.Nevertheless, the crossing swords scheme seems to provide more efficient triangulations even for this generalized problem.Crossing swords as well as generalized K1 is the only schemes which are able to fulfill the given accuracies exactly.However, the number of simplices in generalized K1-triangulations growths much faster than in crossing swords triangulations with increasing accuracy.For ε = 0.05, K1 requires even factor twice as many simplices as crossing swords.The poor adaptability to the prescribed approximation error exhibited by the refinement methods, i.e., red-refinement, longest-edge bisection and max-error refinement, can be seen from the fact that in order to satisfy an accuracy of ε = 0.5, they produce triangulations that have actual approximation errors of 0.1875 and 0.2344, respectively.To visualize the observed results, the triangulations of the considered triangulation schemes for an accuracy of ε = 0.5 are shown in Fig. 10a.Besides the mere number of simplices, it is especially noticeable that the maximum error scheme produces relatively skinny simplices, which can lead to numerical problems in the evaluation of the approximation for small errors.

Conclusion
In this work, we have introduced a novel triangulation scheme called crossing swords to interpolate bivariate products x y over rectangular domains.We showed that a crossing swords triangulation requires at most √ 5/2 times as many simplices as an ε-optimal triangulation for any approximation accuracy ε > 0. Crossing swords thus outperforms all previously known triangulation schemes in the literature, which we have proved theoretically and also underlined by exemplary numerical results.We have also proved that crossing swords triangulations are ε i -optimal under the condition that each simplex has one axis-parallel edge and ε i = 1 /16i with i ∈ N holds.Supported by computational tests, we conjecture that this also holds if the simplices are not required to be axis-parallel.Future research could address the extension of interpolations to general approximations.In [2], it was shown that optimal triangulations of R 2 consist of slightly fewer triangles if, instead of interpolations, approximations with a deviation by a constant factor are considered.A similar result is presumably achievable for box domains.In addition, triangulations for polytopes other than the rectangular box domains are canonical candidates for further investigations.

3 . 1 Lemma 3 . 3 (
to show that we can reduce the triangulation problem over D to the unit box U := [0, 1] × [0, 1] by scaling the prescribed accuracy with the area of the box.In particular, we show how to transform any ε-triangulation T U of U into an νε-triangulation T D of D such that |T U | = |T D | holds, with ν:=( x − x)( ȳ − ȳ) being the area of D. The triangulation T D is obtained by a linear mapping of the nodes N (T U ). Invariance of triangulations under scaling and shifting) Let f U : U → R be a pwl.ε-interpolation of F defined by a triangulation T U of the unit box U .Further, let L : R 2

|T | = 2 :
Relations (2)-(6) reduce the set of possible triangulations to those where |N (T )| = 4, K = 4 and E = 5.Therefore, only two triangulations of U are possible, namely the ones using one of the diagonals to triangulate U .Either diagonal leads to an optimal triangulation as they produce the same approximation error of ε * 2 = 1 /4, attained at their center.We refer to these equivalent triangulations by T * 2 , and one example is shown in Fig. 3a.Next we consider |T | = 3: Relations (2)-(6) only allow for triangulations with |N (T )| ∈ {4, 5}.Assuming |N (T )| = 4, we obtain K = 3, which is a contradiction to Eq. (4).Assuming |N (T )| = 5 instead, we obtain K = 5.The latter implies that all nodes have to be on the boundary of U .Due to symmetry, it does not matter if the additional node is on a facet which is either parallel to the x-axis or the y-axis.

Fig. 3
Fig. 3 Triangulations with |T | ≤ 5; (*) a, b, c and i are optimal triangulations with respect to the number of simplices used

Fig. 9 Corollary 3 . 22 (
Fig. 9 Comparison of triangulation scheme in terms of the number of simplices for prescribed approximation accuracies

Definition 2.5
Consider a triangulation T of a polytope P ⊂ R d and let g : P → R be a pwl.interpolation of a function G : P → R w.r.t. to T .We call