There are various notions about being at the boundary of the feasible set. The most general case is to say that a polytope p is at the boundary of q, if p ⊆ ∂ q. Now, in our case we want a subset to be in the relative boundary of the feasible set. To make it clear, we introduce the notion of border as follows.

Figure 1 illustrates when facets are border or not. We extend the concept of border status of a face defining its border level.

The practical impact of certifying the border level of polytope subsets is given in the following results, which follow directly from the definitions.

Notice that in order to check the border status of a face f using the above definition and lemmas, is complicated. It requires the set of faces φ of q stored with the dimension and a check whether one of them contains face f. For an n-box, the faces are easily defined by lower and upper bounds of the variables. For a simplicial set, the faces are defined by the so-called face graph. Let I : = { 1 , 2 , … , | W | } be the index set of vertex set W, then any subset J ⊂ I corresponds to a face of a simplicial set conv ( W ).

For a non-simplex polytope, this is not the case. It is possible to verify a subset W J : = { w j ∣ j ∈ J } of W is a face φ by the following result. Let c be the centroid of W J . Consider the following LP problem: (4) min ∑ j ∈ J λ j s.t. ∑ j ∈ I λ j w j = c , ∑ j ∈ I λ j = 1 , λ j ⩾ 0 , j ∈ I .

Using Definition 2 directly would imply using (4) and storing all faces φ of q with their dimension. Checking a face f ⊂ φ requires solving (4) taking the centroid c of f and the subset J which defines φ. Checking this for all faces φ does not sound practical. There are some simple rules that can be derived from properties.

As discussed in Hendrix et al. (2024), there are several ways to refine the partition, i.e. to split the subsets further. In this paper, we will focus on bisecting the subset along the widest coordinate of its interval hull, following traditional methods in iBB. To avoid dealing with the face set, we will keep track of border levels with a focus on facets and edges of a subset. Notice that all faces of the feasible set q are border. The first polytope subset is p = q. Its facets f have border level dim ( f ) = dim ( q ) − 1. In our procedure, we will keep track of the border status of facets and edges of polytope subsets p in the division (see Section 4.2), and update them when p is reduced, see Section 4.6.

Consider two points x , y ∈ p, x ≠ y that define a direction from x: d = y − x, over polytope p.

As we consider continuously differentiable functions, monotonically increasing follows from d T ∇ h ( x ) > 0, Hendrix and G.-Tóth (2010). Let ∇ h ( p ) be the range of ∇ h ( x ), ∀ x ∈ p. Then d T ∇ h ( p ) : = [ d T ∇ h ( p ) _ , d T ∇ h ( p ) ‾ ] is an interval with bounds containing all directional derivative values d T ∇ h ( x ), x ∈ p where d is a direction in p. In this study, we use Interval Automatic Differentiation over □ p, to get bounds of the derivatives of the objective function over p, ∇ h ( p ) ⊆ g ( □ p ) : = [ g _ , g ‾ ], see Moore et al. (2009), Rall (1981). Then, the corresponding directional derivative d T ∇ h ( p ) is in the range of the interval inner product (5) d T ∇ h ( p ) ⊆ d T g = [ d T g _ , d T g ‾ ] = [ ∑ i = 1 n min { d i g _ i , d i g ‾ i } , ∑ i = 1 n max { d i g _ i , d i g ‾ i } ] . Notice that overestimation (a less tight enclosure) of p by □ p is larger as p deviates more from □ p. This increases the overestimation of d T ∇ h ( p ) in d T g. Let d d : = d T g denote the interval directional derivative enclosure. If 0 ∉ d d, then d is a monotonic direction on p.

In Casado et al. (2025), we focus on the question how to find monotonic directions in sBB. The most important impact is that, as for simplices (Casado et al., 2022), the relative interior of the full or non-full dimensional polytopal partition set cannot contain an optimum.

The main consequence of Proposition 3 is that a non-interior subset p can be replaced by the collection of its facets, which basically have a lower polytope dimension and fewer vertices. This is mainly a theoretical observation. Keeping track of the facets of a polytope is complex. We will come back to this issue. One typical observation following Corollary 1 is the following.

We now focus on the direction from a vertex of p pointing towards the relative interior, d = y − v, v ∈ V.

The argument is that the relative boundary of f consisting of faces not containing v may contain a global minimum point. However, the global minimum point is enclosed in the relative boundary of other facets that do not include v. One should take care in not applying this concept for several vertices simultaneously in all cases.

Figure 2 illustrates the consequences of Corollary 4.

There are several ways to find monotonic directions. In Hendrix et al. (2024), we extended the earlier analysis of Casado et al. (2021) for simplicial partition sets to polytopal partition sets. An easy way to start is to focus on the directions d = c − v from a vertex to the centroid and the directions between vertices. In Hendrix et al. (2024), we show that the best direction in sense of highest or lowest directional derivative bound can be found solving a linear program. This may not be computationally efficient as discussed in Casado et al. (2025).

In this study, we describe a polytope p = { V , E , F , m } as a set of vertices, edges, facets and its dimension, respectively. By performing a bisection, we generate two polytopes, one at left p ℓ = { V ℓ , E ℓ , F ℓ , m }, and one at right p r = { V r , E r , F r , m } of the midpoint of the widest coordinate of □ p. One of the tasks to perform is to determine the cutting facet f c . Determining the vertices of f c is not difficult, but to determine its edges can be a challenge, because a pair of vertices in f c might not be connected by an edge. We first show the difficult case and later on the easier ones in determining f c .

As illustration, Fig. 3 shows some cases of the cutting facet f c in red for a cube bisection. The cubes may be rotated to illustrate a component-wise division where internal division is not in the middle of the widest component. Notice that in general we divide edges, but there are cases where we do not, see left-down graph of Fig. 3.

Table 2 shows the notation for the used elements and the sets to be determined as starting step before performing a bisection. Facets of a facet f j of a polytope p have to be determined only if f j is treated as a polytope, which will occur only after a reduction to it, see Section 4.5. Therefore, we do not need to identify the facets of a facet in the division of p, but only the set of vertices and edges a facet has.

The difficulty is illustrated by the labelling provided in Table 2. After the identification of the sets, the aim is to divide edges at E d generating new vertices at c that are added to V c . For each divided edge, we generate a new vertex at c and two new edges that are stored at new sets: •

E d ℓ : set of edges at left of c, as left result of division of edges in E d .

In the case where the facet to be divided has no edges to divide (see top and bottom facets of the cube at left-down graph of Fig. 3), we have to separate edges generating two facets; one with edges at left and store them in F d ℓ and one with edges at right and store them in F d r . At this point, we have open facets sets at the left F d ℓ and right F d r . Figure 4 shows an example of sets F d ℓ and F d r of a cube bisection.

What remains to identify are edges E c of f c after construction of the complete vertex set V c . The difficult part of a bisection algorithm based on only vertices is that it requires checking, whether a pair of vertices ( u , v ) is an edge of a polytope, in our case for f c . In López and Duval (2012), one can find algorithms to compute edges from a vertex set. They are based on finding supporting hyperplanes using linear programming (LP). A simple procedure is to check by LP (4) whether the midpoint c = v ℓ + v k 2 of two vertices v ℓ and v k is on a facet of p. Such procedures require running an LP routine. This does not seem attractive in an iterative process for each pair of vertices of f c in a spatial B&B search. Moreover, in our experience, LP routines are getting sensitive to numerical errors when the subsets are getting small.

We use a different approach in Algorithm 2 to determine the edges of k-facet f c , for the most complicated cases, | V c | > 3, k > 1. Easier cases will be described later in this section, see graphs of Fig. 5. Algorithm 2 uses a set F d r to determine all edges e ∈ E c . This could be F d ℓ as well, because both need the same edges of facet f c . It starts to determine the set of vertices at c for each facet f j ∈ F d r . We also need to keep track of the indices of the facets the vertices are in. Thus, we use a set R with elements r = ( U , I ) ∈ R having two components: U, the set of vertices at c and I, the set of indices of the facets the vertices at U belong to. In line 2 of Algorithm 2, we store the vertices at c and the index of the facet for each f j d r ∈ F d r in R. The main computational burden of the procedure is in the while in line 6 of Algorithm 2. This computational burden increases with the dimension of f c . It is based on the iterative intersection of vertices U and the union of the corresponding indices of the facets I for elements in R, while at least two vertices are in U for each r ∈ R. Intersections with less than two vertices are not considered (they can not generate an edge) and we also have to take care to avoid duplicates in the process. This process can be seen as the traversal of the face-graph of f c , having f c at the top level, and its facets (determined by their vertices) in the next lower level. The process moves down, level by level, until the level with edges is reached. Notice that the while loop will not be executed in the illustrative example in Fig. 4, because each U j has exactly two elements, so we find all edges right away. Finally, Algorithm 2 adds new edges to the corresponding open facets at F d ℓ and F d r to close them. Then, F d ℓ is added to F ℓ ⊂ p ℓ and F d r to F r ⊂ p r . Facet f c completes both F ℓ and F r . The bottom graphs of Fig. 4 show an example of resulting sets p ℓ and p r .

As mentioned above, there exist easier cases to determine cutting facet f c . If the set of facets to be divided F d is empty (see right graph of Fig. 5), vertices in V c and edges in E c of f c are identified. In this case, the set of facets in p ℓ is F ℓ ∪ f c and the set of facets in p r is F r ∪ f c . When the number of vertices in V c after dividing edges is three, they must be all connected in f c , see left graph of Fig. 5. So, Algorithm 2 is not needed.

For 2-polytopes, it is easy to determine cutting facet f c , because f c is a segment connecting two vertices, which can be new or not. Figure 6 shows some examples.

Finally, if p is 1-polytope (segment) the cutting 0-facet is the new vertex.

For keeping track of the border level, we can exploit the following rules: •

Facets f and edges e of feasible set q, an m-polytope, are border and have border-level b l ( f ) = m − 1 and b l ( e ) = 1.

Interval Arithmetic and Automatic Differentiation are easy choices to get bounds of the objective and its derivatives over a box. In a very summarized way, Interval Arithmetic replaces reals by intervals and real operations by interval ones. The very first objective of Interval Arithmetic was to avoid rounding errors of Computational Arithmetic. For polytopes (and simplices) we can use the subset interval hull to get bounds of the objective function over p: h ( p ) ⊆ h ( □ p ) and its derivatives: ∇ h ( p ) ⊆ g ( □ p ). So, iBB tools can be used in pBB.

The so-called centred form of the Taylor extension may improve the bounds of the direct extension from real arithmetic to intervals h ( □ p ) by using first-order derivatives h c ( □ p ) = h ( c ) + g ( □ p ) ( □ p − c ), for some c ∈ □ p, see G.-Tóth et al. (2021) for more details on bounding rules. We will use the evaluation of the centre of p in h c ( □ p ), instead of the centre of □ p, which usually obtains a smaller maximum distance from centre to vertices, which is required in the h c ( □ p ) computation.

As a selection rule, we choose the polytope subset p ∈ Λ with the lowest updated lower bound h _ ( □ p ) = max { h _ ( □ p ) , h c _ ( □ p ) } to be processed in the next iteration of Algorithm 1.

Using the centred form, the centre of p is evaluated. It can update the best upper bound of the solution so far, the so-called incumbent represented by h ‾. Symbol h _ denotes the lower bound of the solution, which is h _ ( □ p ) for the selected polytope.

Algorithm 1 terminates when h ‾ − h _ ⩽ α. The interval [ h _ , h ‾ ] and the remaining inclusion set Λ, which contains the minima, are the output of Algorithm 1.

The theoretical results give rise to concrete tests that can be used to reject subsets from further consideration or to reduce them to lower dimensional subsets. First of all, the RangeUp test rejects a subset p if h _ ( □ p ) > h ‾, see Section 4.4. If the test fails, we can evaluate the centred form to test h c _ ( □ p ) > h ‾. The application of the RangeUp test to already evaluated subsets at Λ is the so-called CutOff test.

We now focus on exploiting the theoretical monotonicity properties discussed in Section 3.2 to decide on rejection and reduction. A necessary condition to reject or reduce a polytope p based on monotonicity is 0 ∉ g ( □ p ). It is important to differentiate if the m-polytope p is full dimensional ( m = n) or not ( m < n), and if it has border facets or not. For m-polytopes, m < n, condition 0 ∉ g ( □ p ) is not sufficient and we need to find a directional derivative 0 ∉ d d = d T g ( □ p ) in p, see (5). We will use the notation d d + for d d _ > 0 as positive directional derivative and d d − for d d ‾ < 0. In this study, we limit the search to directions from vertices of p, d d v , v ∈ p to the centre and to other vertices of p. Figure 7 shows a hypothetical example where both d d v − and d d v + exist in v. We are interested mainly in d d v − , see Corollary 4. Hence, negative directional derivative d d v − could replace the information about monotonic directional derivative in v when it is d d v + . Notice that a positive directional derivative d d v + from v to vertex u implies having a negative d d u − .

The rejection and reduction cases when 0 ∉ g ( □ p ) are as follows.

In case 2(b)iiB, we can stop after one negative direction d d v − is found or look for all of them and apply Corollary 4 using the vertex that removes more border facets. If all border facets are removed, p is not reduced but rejected. If p is a 1-polytope (segment), it can only be reduced to a vertex of the initial feasible set q, which is a border 0-facet of p.

Consider reduction of subset m-polytope p = { V , E , F , m } to a border facet f j ∈ F. Then, f j is evaluated, see Section 4.3. Before storing f j in Λ, the RangeUp test is checked. If it passes the test, we have to store f j as a polytope, f j = { V j , E j , F j , m − 1 }. This requires to determine the set of ( m − 2 )-facets f i f j ∈ F j of f j . Each f i f j is obtained by a nonempty intersection of polytopes f j ∩ f i , j ≠ i, f j , f i ∈ F ∈ p.

We observed in the experiments that when selecting the next reduced-evaluated-RangeUp p from Λ, it is worth to try to reduce p again with the new derivative bounds, instead of directly dividing it. We do not apply more than one reduction-evaluation-RangeUp before to store p in Λ, because this would affect the selection rule (see Section 4.4), which might increase the computational burden of the algorithm. However, by storing p after a reduction-evaluation-RangeUp, the CutOff could remove p from Λ in a next iteration.

The determination of the border level of facets and edges after a polytope division was discussed in Section 4.2. Lemma 1 determines the border status (border level) of a facet. Facets generated by successive division of a facet f maintain the border status of f. All facets of feasible set q are border. A first non-border cutting facet f c is generated by the bisection of the widest component of □ q. The algorithm will never reduce to a cutting facet, or facets generated by a successive division of a cutting facet, because they are not border and are labelled as non-border after a division.

Consider an m-polytope p = { V , E , F , m } to be reduced to a border ( m − 1 )-facet f j = { V j , E j , F j , m − 1 }. As Section 4.5 discussed, each ( m − 2 )-facet f i f j of f j is obtained by the non empty intersection f j ∩ f i , j ≠ i, f j , f i ⊂ F ⊂ p. To determine the border status of the facets f i f j of f j , Proposition 2 states that if both f i and f j are border, then f i f j is border.

The next question concerns the border status of the facets f i f j = f j ∩ f i after reduction of p to f j and f i is not border due to being a cutting facet or a division of it. The problem is that face f i f j ⊂ p is a facet of both f j and f i , and it can be border. Examples are provided at the bottom left graph of Fig. 3 and right graph of Fig. 5, where there exist border edges in f c . In general, f i can have a border face. So we need some tools to determine its border status. Proposition 5 characterizes a non-border facet without border faces.

Proposition 5 also shows that facet f c under this conditions does not contain a face of p. If E c = ∅ in the input of Algorithm 2 and f c has no vertex of q, then f c has no border faces. Therefore, f i f j and any of its subsets generated by division are not border, as well.

However, if f c contains faces of p, the problem is to identify them and to determine they are not border and if they are also in another border facet. In our algorithm, we try to simplify this task using the condition of Proposition 6.

Proposition 6 is only relevant for a face f i f j after a reduction, because all facets were already labelled in a division. It provides a necessary, but not sufficient test for facet f i f j to be border. If the condition is not fulfilled, the facet is not border.

Moreover, Proposition 2 says that f i f j is border if both f i and f j are border. In the algorithm, we test the condition of Proposition 6 on f i f j only when f i f j was not labelled as border by Proposition 2 and it could not be labelled as non-border by Corollary 5.

Proposition 6 requires to keep track of the edge border levels after reducing to a border facet f j . Rules in Section 4.2 show that the border level of an edge can be inherited from a divided edge or from the facet it is generated in. By reducing an m-polytope p to a border ( m − 1 )-facet f j , we diminish the dimension by one. This implies that, when m > 3, we also have to decrease the border level of edges of p with b l ( e ) ⩾ ( m − 1 ), before they are inherited by the border facet f j . Figure 8 shows an illustrative example.

Figure 8 provides an example of Proposition 5 in the red cube f c . Figure 8 also aids to observe that Proposition 6 provides a necessary, but not sufficient condition to label a facet f i f j of the border facet f j to which we reduce. All 2-faces of f c are not border. They are in the relative interior of the corresponding 3-facet of q. Consider the 4-cube q ℓ in the middle graph. It has eight 3-cube facets. One of them is the cutting facet of q in red. If we reduce q ℓ to the 3-cube facet f j at the left, we have to determine its facets f i f j by intersecting f j with the other 3-cube facets of q ℓ . All facets of q ℓ are border apart from the red 3-cube facet f i . So, all but one f i f j are border and we have to determine the border status of f i f j = f j ∩ f i , which is the red square in the right graph. Face f i f j is in the relative interior of q and satisfies the condition of Corollary 5. So, it is labelled as non-border in pBB because f c has no border faces. Although, Proposition 6 will not be checked by pBB, it is satisfied, i.e. there exists one vertex in two edges with border level 2. So, Proposition 6 does not provide a test to guarantee the facet is border.

In case the condition of Proposition 6 is not satisfied, we label f i f j as border. The open question is whether a non-border f i f j exists not satisfying Corollary 5 and satisfying Proposition 6. In such case, we are labelling f i f j as border incorrectly and the algorithm may reduce to it. However, pBB still converges to the global minimum, because reduction is only to border facets and facets that might contain a global minimum are not discarded in the search. But pBB may evaluate polytopes that do not need to be tested.

Each algorithm has a different input for the same polytope q as feasible set:

In order to compare the three spatial B&B algorithms as fair as possible, the following B&B rules were taken for all algorithms.

The feasible area for all algorithms is the same. Some considerations about feasibility of a box in iBBLC are the following. A point is treated as a degenerated box in its feasibility evaluation. Having linear constraints A x ⩽ b, for each slack s j = A j x − b j , over a box x, we have:

The linear constraints for the instances are provided in Appendix A. For the iBBLC algorithm we show if the minimum point x ∗ is feasible (F), infeasible (I) or undetermined (U). Additionally, we show the exponent of the distance of s j to zero in scientific notation. Examples of this distance (which is different from interval distance to a point) are shown next. For instance, e-x in Ie-x is the exponent of the very first infeasible s j given in s j _. In Fe-x, e-x is the minimum exponent for all feasible s j , given in some s j ‾. If there exists one s j ‾ = 0.0, we show F0 instead. In Ue-x, e-x is the minimum exponent for all limits s j _, s j ‾ of undetermined (might exist feasible s j ) slacks s j .

Due to computational representation of numbers, two numbers are considered equal if their absolute difference is less than 10 − 12 . This is used to determine if a new vertex (point) already exists and has already been evaluated. This is also relevant for polytope division in order to determine if a cutting vertex already exits. Special care is required to the storage of a set S in Λ in an efficient way. We use an AVL tree, sorted by h ( S ). A linked list is used in an AVL tree node for equal values of different h ( S ). In the list, we use LIFO (last in, first out) to be more efficient in insertion and extraction. Interval x < y when x _ < y _ or ( x _ = y _ and x ‾ < y ‾). Vertices are also stored in an AVL tree sorted by their coordinates. Thus, no duplicates exist. We can remove a vertex when it is no longer in a set (to save memory) or let it stay forever (fewer vertex evaluations but more memory requirement). For this experimentation, we have chosen not to remove them.

The considered instances are described in Appendix A. For sBB, non-simplicial polytopes are partitioned into simplicial subsets using the algorithm in Assarf et al. (2017). To use sBB for more than one simplex in the initial list of subsets, we had to extend it to deal with several simplices where not all facets are border initially. One should take care because, in such a case, sBB may discard a non-border facet which includes the solution, as the algorithm might find a negative d d v − to that non-border facet, see G.-Tóth et al. (2024).

Table 4 shows the results of the three algorithms for the designed set of instances. For non full dimensional instances 2Pol and 4Pol, see Appendix A.1 and A.3, no box generated by the iBBLC algorithm is feasible and the number of boxes increases significantly with the dimension. Although the execution time for all algorithms is less than 0.01s, iBBLC generates 150 times more subsets than sBB and pBB for the 2Pol instance. The differences are larger for the 4Pol instance: iBBLC could not solve it in less than 5 minutes and sBB generates thousands of simplices, but pBB only hundreds of polytopes.

Table 5 describes the meaning of the colours used in the illustrative Figures 9 to 13.

Figure 9 shows graphically the iBBLC execution over instance 2Pol. Yellow lines represent the linear constraints. There exist boxes rejected by infeasibility that are graphically overlapping the constraint, but they are actually not doing so. This visual effect is due to the graph only showing x 1 , x 2 coordinates and there is more than one box sharing x 1 and x 2 components, see Fig. 14 in Appendix A.1. Even when no boxes in iBBLC are feasible, at least some evaluated centres of boxes, by the interval centred form, are feasible. Finding a feasible point is very important in constrained problems in order to use the RangeUp and CutOff tests. The refinement generates a lot of infeasible and undetermined boxes before reaching the termination criterion.

Figure 10 shows the steps of the sBB algorithm on the simplicial partition of the 2Pol instance. Consider the graphs in a row-wise order, where in each graph the simplex and its interval hull are drawn. In the first graph, the simplex is rejected by a negative directional derivative d d v − to a non-border facet. In the second graph, the next simplex is reduced to the left vertex, because a negative directional derivative d d v − was found for the other two vertices. In the third graph, a simplex is reduced to the border edge on the left, because there exists a d d v − from the vertex not belonging to the border edge. In the fourth graph, the reduced edge is evaluated and stored. In the fifth graph (second row), the next selected simplex is reduced to the vertex solution v ∗ , the leftmost vertex. Notice the d d v ∗ + is negative in the other direction. In the sixth and seventh graphs, selected and evaluated simplices were rejected by the RangeUp test. In the last graph, the stored border edge is selected and reduced to v ∗ which was already evaluated. Reduced and evaluated simplices are not directly divided, but checked for additional reduction first. This reduces the number of simplex evaluations.

Although we use □ s to get bounds of the objective over the simplex, which incurs in overestimation of the volume of the simplex, the search is restricted to the feasible area. This avoids the feasibility test as is done in iBBLC and its accuracy issue.

Figure 11 illustrates the application of algorithm pBB to the 2Pol instance. Consider the vertices clock-wise, starting at v 1 at the rightmost vertex. The left graph corresponds to two negative directional derivatives d d v − in red from vertices v 1 and v 2 to the centre. This implies removing the interior of 2Pol and the facets containing these vertices. We can use more than one vertex, because we are in a 2-polytope, see Fig. 2. In the middle graph, vertices v 1 , v 7 and v 8 have a negative vertex to vertex direction d d v − shown in red. Direction v 5 ∗ to v 4 results in d d v + , which is negative in the opposite direction. After checking all directional derivatives, the initial polytope is reduced to edge ( v 5 , v 6 ), see Section 4.5. The right graph illustrates the reduction of that edge to the solution vertex v 5 ∗ . For this instance, pBB does not perform a division, only reductions. Therefore, the number of evaluated polytopes is just two and the number of evaluated points is just three.

For the full dimensional 3Pol case in Table 4, we use five different objective functions generating five instances, see Appendix A.2 and Fig. 15 in Appendix A, given by: