Out-of-sample estimation for a branch-and-bound algorithm with growing datasets

Sass, Susanne; Mitsos, Alexander; Nikolov, Nikolay I.; Tsoukalas, Angelos Article — Published Version Out-of-sample estimation for a branch-and-bound algorithm with growing datasets Journal of Global Optimization Provided in Cooperation with: Springer Nature Suggested Citation: Sass, Susanne; Mitsos, Alexander; Nikolov, Nikolay I.; Tsoukalas, Angelos (2025) : Out-of-sample estimation for a branch-and-bound algorithm with growing datasets, Journal of Global Optimization, ISSN 1573-2916, Springer US, New York, NY, Vol. 92, Iss. 3, pp. 615-642, https://doi.org/10.1007/s10898-025-01514-4 This Version is available at: https://hdl.handle.net/10419/323681 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. http://creativecommons.org/licenses/by/4.0/ Journal of Global Optimization (2025) 92:615–642 https://doi.org/10.1007/s10898-025-01514-4 Out-of-sample estimation for a branch-and-bound algorithm with growing datasets Susanne Sass1·Alexander Mitsos1,2,3 ·Nikolay I. Nikolov4,5 · Angelos Tsoukalas6 Received: 10 May 2024 / Accepted: 12 June 2025 / Published online: 23 June 2025 © The Author(s) 2025 Abstract In [Sass et al., Eur. J. Oper. Res., 316 (1): 36 – 45, 2024], we proposed a branch-and-bound (B&B) algorithm with growing datasets for the deterministic global optimization of parameter estimation problems based on large datasets. Therein, we start the B&B algorithm with a reduced dataset and augment it until reaching the full dataset upon convergence. However, convergence may be slowed down by a gap between the lower bounds of the reduced and the original problem, in particular for noisy measurement data. Thus, we propose the use of outof-sample estimation for improving the lower bounds calculated with reduced datasets. Based on this, we extend the deterministic approach and propose two heuristic approaches. The computational performance of all approaches is compared with the standard B&B algorithm as a benchmark based on real-world estimation problems from process systems engineering, biochemistry, and machine learning covering datasets with and without measurement noise. Our results indicate that the heuristic approaches can improve the final lower bounds on the optimal objective value without cutting off the global solution. Aside from this, we prove that resampling can decrease the variance of the lower bounds calculated based on random initial datasets. In our case study, resampling hardly affects the performance of the approaches which indicates that the B&B algorithm with growing datasets does not suffer from large variances. Keywords Nonlinear programming ·Spatial branch and bound ·Parameter estimation · Overfitting ·Resampling BAngelos Tsoukalas [email protected] 1Process Systems Engineering (AVT.SVT), RWTH Aachen University, 52074 Aachen, Germany 2JARA-CSD, 52056 Aachen, Germany 3Institute of Climate and Energy Systems: Energy Systems Engineering (ICE-1), Forschungszentrum Jülich GmbH, 52425 Jülich, Germany 4Institute of Statistics, RWTH Aachen University, 52056 Aachen, Germany 5Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria 6Department of Technology and Operations Management, RSM Erasmus University Rotterdam, 3062 PA Rotterdam, Netherlands 123 616 Journal of Global Optimization (2025) 92:615–642 1 Introduction The validity of parameter estimation results increases significantly when avoiding both suboptimal solutions and overfitting. Suboptimal solutions can be excluded when using deterministic global optimization (DGO) methods [1,2]. Common DGO methods [3,4] like the branch-and-bound (B&B) algorithm [5,6] reach their limitations when solving large-scale nonconvex optimization problems arising from many parameter estimation problems [7–9]. To overcome these limitations, we [10] extended the standard B&B algorithm by using growing datasets. Briefly speaking, we start the B&B algorithm with a reduced model obtained by picking only a small subset of the full dataset and augment the dataset until converging to the full dataset provided. For this, we introduced augmentation rules which decide for each processed node whether to augment the dataset or to branch the parameter domain. We have proven that the B&B algorithm with growing datasets remains a DGO method, see Theorem 2 of [10]. A key point of the proof is that the lower bounding problem based on reduced datasets yields a valid lower bound for the original problem, i.e., the parameter estimation problem based on the full dataset. Another key point is the choice of an augmentation rule which guarantees reaching the full dataset eventually. However, when data reduction allows for a much looser lower bound and objective value, convergence may be slowed down or even prevented depending on the augmentation rule, see Example 1 in [10]. This brings us back to the topic of overfitting. In the context of machine learning, overfitting refers to the phenomenon of models performing well on the data used for training but fail to generalize to new data [11,12]. In practice, the actual performance of the trained model is therefore commonly estimated based on independent datasets, the so-called validation sets. We transfer these findings to the B&B algorithm with growing datasets by using out-of-sample evaluations for improving the lower bounds calculated based on the reduced dataset (reduced lower bound). In detail, we use an approximated lower bound given by a combination of the reduced lower bound and the out-of-sample evaluation for (i)a novel augmentation rule and (ii)for pruning. While (i)yields a heuristic rule aiming to improve the performance of the deterministic approach, with (ii)we loose the theoretical guarantee for converging towards the global solution and therefore obtain mere heuristic approaches. To estimate whether the heuristic pruning indeed cuts off the global solution, we propose a post-processing check of the final lower bound after the termination of the B&B algorithm. In [10], we proposed to pick the data points for the reduced datasets randomly from the full dataset. In this article, we exploit this randomness for improving the reduced lower bound for both the heuristic augmentation rule and the heuristic pruning. In fact, resampling is a common tool in statistics to reduce the bias and variance of estimators, cf. [13,14]and[15, Chapter 5]. In our approach, resampling a reduced dataset comes with the calculation and optimization of a lower bounding problem, making it computationally costly. We therefore propose to use two subsamples for updating the reduced lower bound based on the initial dataset. All proposed extensions of the B&B algorithm with growing datasets are implemented in our open-source solver MAiNGO1[16] which is a DGO solver for factorable mixedinteger nonlinear programs. We perform an extensive case study comparing the deterministic and heuristic approaches of the B&B algorithm with growing datasets with the standard B&B algorithm as a benchmark. For this, we revisit real-world applications from our previous work and collected further models from literature resulting in 13 different parameter 1Available at https://git.rwth-aachen.de/avt-svt/public/maingo. 123 Journal of Global Optimization (2025) 92:615–642 617 estimation problems covering the fields of process systems engineering, biochemistry, and machine learning. Ready-to-use implementations of the problem formulations are published in our open-source repository GloPSE2. The remainder of this article is structured as follows. Section 2comprises the mathematical background, where we start with an overview on our notation in Section 2.1. Subsequently, the main algorithmic and theoretical advances of the B&B algorithm with growing datasets from [10] are recalled in Section 2.2. In Section 3, we obtain an approximate lower bound for the original problem based on the reduced lower bound and out-of-sample evaluation. This approximation is used for a novel augmentation rule in Section 3.1, and for introducing a heuristic approach including a post-processing step in Section 3.2. In Section 3.3, we transfer the findings to estimation problems minimizing the mean squared error yielding a second heuristic approach. In Section 4, we propose a heuristic lower bound based on resampling. In Section 5, we perform an extensive case study evaluating the proposed post-processing procedure and comparing all approaches of the B&B algorithm with growing datasets with the standard B&B algorithm, before we conclude in Section 6. 2 Preliminaries 2.1 Problem formulation and notation As in our previous work [10], we focus on finding globally optimum parameter values p∈P minimizing min p∈P (xd,yd)∈D g(p;xd,yd)(SSE) s.t. h(xd,yd;p)≤0∀(xd,yd)∈D  h(p)≤0,(PE) where g(·; xd,yd)=(f(xd;·)−yd)2is the squared prediction error of a parameterized model function f(·; p):Rm→Rfor data point (xd,yd)∈DR m×R,andh(xd,yd;·) and  hare the residual functions of the data-dependent and data-independent inequality constraints on parameters p∈P, respectively. Let the parameter domain PR nbe a closed, bounded box. In case of integer parameters, let the parameter domain consist of subsequent discrete values. Note that g(·; xd,yd)is nonnegative, and strictly positive in case of model-data mismatch caused by measurement errors, suboptimal parameter values, or model misspecification. In the following, we call Dthe full dataset and parameter estimation problem (PE)the original problem.Letgcv(·; xd,yd):P→Rbe a non-negative convex underestimator of g(·; xd,yd)as well as hcv(xd,yd;·)and hcv be convex underestimators of h(xd,yd;·)and h, respectively. Then, the optimal solution of the convex optimization problem min p∈P (xd,yd)∈D gcv(p;xd,yd)(LBP) s.t. hcv(xd,yd;p)≤0∀(xd,yd)∈D  hcv(p)≤0, 2Available at https://git.rwth-aachen.de/avt-svt/public/glopse. 123 618 Journal of Global Optimization (2025) 92:615–642 gives a lower bound on the optimal solution of the original problem, which we will call full lower bound.Werefertoreduced problem and reduced lower bound when replacing the dataset with a reduced dataset Dk⊆Dwithin (PE)and(LBP), respectively. Note that, in common DGO software, the lower bounding problem (LBP) is further relaxed to a linear program without losing its validity for the original problem by linearizing valid convex underestimators. In the B&B node Nk=(Pk,Dk)processed in B&B iteration k∈N, optimization problems (PE)and(LBP) are solved for parameter domain Pk⊆Pand dataset Dk⊆D.Letuk be the best upper bound found until iteration k.Letlfull kand lred kbe the full and reduced lower bound calculated in node Nkwith optimal parameter values plb,full kand plb,red k, respectively. If a lower bound calculated based on a reduced dataset cannot be guaranteed to be valid for the original problem, we will call it heuristic lower bound lkcompared to valid lower bounds lk. 2.2 The B&B algorithm with growing datasets As proposed in [10], the B&B algorithm with growing datasets extends the standard algorithm by a data reduction step before entering the B&B loop and a subroutine for deciding whether to branch or augment a node, see Figure 1. For this, we associate all nodes Nk=(Pk,Dk)∈N of the B&B tree with both a parameter domain PkPand a dataset Dk⊆D.Inthe subroutine, cf. Subroutine 1 of [10], we call an augmentation rule A:N→{True,False}. If A(Nk)=True, the dataset is augmented, i.e., we add one child node Nnew =(Pk,Dnew) with Dnew Dk. Otherwise, we branch the parameter domain while retaining the dataset, i.e., we add two child nodes Nnew,1 =(Pk,1,Dk)and Nnew,2 =(Pk,2,Dk)with partition Pk,1∪Pk,2=Pk. With this, we obtain an algorithm which is guaranteed to converge towards the global optimum of (PE) if we use a finitely convergent augmentation rule [10]. Definition 1 (Definition 1 of [10]) Let Dk⊆Dbe the dataset used in node Nk,whichis processed in iteration kof the B&B algorithm with growing datasets depicted in Figure 1. The augmentation rule A:N→{True,False}completes finitely, if for any infinite nested sequence of nodes {Nkj}j→∞ with Nkj=(Pkj,Dkj)it holds ∃J<∞:Dkj=D∀j≥J. For completeness of contents, we repeat the main theoretical results of [10]. The interested reader may refer to the original publication [10] for more details and the proves. The theoretical results are based on the following assumptions. Assumption 1 (Assumptions 1 and 2 of [10]) (i) Let f(xd;·):P→Rand h(xd,yd;·):P→Rfor any fixed (xd,yd)∈Das well as  h:P→Rbe continuous. (ii) Let gcv(·; xd,yd)be any nonnegative convex underestimator of g(·; xd,yd)over P for any fixed (xd,yd)∈D. (iii) Let hcv(xd,yd;·)and  hcv be any convex underestimator of h(xd,yd;·)for any fixed (xd,yd)∈Dand  h, respectively, over P. With the help of the nonnegative convex underestimators, cf. Assumption 1(ii), we can construct valid lower bounds based on a reduced dataset. Lemma 1 (Lemma 1 (i) of [10] (adapted)) Let Assumption 1hold. Then, (xd,yd)∈Dred gcv(p;xd,yd)is a convex underestimator of both (xd,yd)∈Dred g(p;xd,yd)and (xd,yd)∈Dg(p;xd,yd)over Pfor any Dred ⊆D. 123 Journal of Global Optimization (2025) 92:615–642 619 Fig. 1 Flow chart of the B&B algorithm with growing datasets [10] where the extensions of the standard B&B algorithm are highlighted with dashed boxes. Finally, we extended the proof of convergence of a standard B&B algorithm given in Theorem 5.26 of Locatelli and Schoen [17]. In particular, we build upon their definitions for exactness in the limit,theisotonic property of underestimators, and exhaustiveness of branching given in Definition 5.4, Equation (5.34), and Definition 5.5, respectively, of [17]. Theorem 1 (Theorem 2 of [10] (adapted)) Let Assumption 1hold. We apply a spatial B&B algorithm with optimality tolerance ε>0to (PE). Let the convex underestimators gcv D,h cv(xd,yd;·)∀(xd,yd)∈D, and  hcv be exact in the limit. Let gcv(·; xd,yd)and 123 620 Journal of Global Optimization (2025) 92:615–642 Table 1 Datasets, corresponding optimal objective and scaled optimal objective as given for Example 1 of [10] Data points Opt. obj. Opt. obj. scaled with |D| |Dred| D0={(0,0.6)}0.0 3 ·0=0 D1={(0,0.6), (0,1)}0.08 3 2·0.08 =0.12 D={(0,0), (0,0.6), (0,1)}0.5067 0.5067 hcv(xd,yd;·)satisfy the isotonic property ∀(xd,yd)∈D. Let the subdivision process of the B&B algorithm be exhaustive. If we use an augmentation rule Awhich completes finitely, then the B&B algorithm with growing datasets depicted in Figure 1terminates after a finite number of iterations and •either establishes that the problem is infeasible if the final upper bound equals infinity •or returns an ε-optimal solution if the final lower bound is finite. Note that it is essential for the proof of finite convergence that the augmentation rule completes finitely. Otherwise there may be a gap between reduced and full lower bound preventing augmenting or even pruning based on the reduced dataset, see the following example. Example 1 (Example 1 of [10] (adapted)) We want to solve min a∈[0,25](a−1)2+(a−0.6)2+(a−0)2, where we find a linear function f(x;a)=a·xthrough the origin and three data points D={(1,0), (1,0.6), (1,1)}at the same input x=1, with a naive implementation of the B&B algorithm with growing datasets using augmentation rule Ascaling with ρ=1. Assume that the (reduced) datasets are chosen as in Table 1. Even if the lower bound calculated based on reduced datasets D0and D1is exact, i.e., equals the respective optimal objective at the optimum point p∗, gaps would remain  (xd,yd)∈D0 gcv(p∗;xd,yd) (xd,yd)∈D1 gcv(p∗;xd,yd)  (xd,yd)∈D gcv(p∗;xd,yd)≤ (xd,yd)∈D g(p∗;xd,yd) Branching does not help either since we already assume exact lower bounds. With that, neither augmenting nor convergence is possible.  3 Out-of-sample estimation with growing datasets The reduced lower bound is valid for the original problem but may be too loose for augmenting or pruning, cf. Section 2.2 and [10]. In this case, the time savings due to the data reduction may be canceled out by the unnecessary large B&B tree. Thus, we propose to apply out-of-sample evaluation for a better approximation of the full lower bound. We expect that the solution of the reduced bounding problem plus this additional evaluation remains computationally faster 123 Journal of Global Optimization (2025) 92:615–642 621 than calculating the full lower bound. In detail, we evaluate (LBP) over all remaining data points D\Dkat the optimal solution point plb,red kyielding the heuristic out-of-sample lower bound  loos k:=  (xd,yd)∈D\Dk gcv(plb,red k;xd,yd). Note that we define the out-of-sample lower bound only for proper subsets DkD. Based on empirical evidence, we expect the out-of-sample lower bound to be valid for the original problem. In fact, the out-of-sample lower bound sums up less nonnegative terms than the full lower bound, where each term is a function evaluation at point plb,red kwhich is typically close or equal to the optimal solution point plb,full k. Postulate 1 For B&B nodes with reduced datasets DkD, out-of-sample lower bound loos k is smaller than or equal the full lower bound lfull kand therefore valid for the original problem. Postulate 1may be violated for pathological cases as shown in the following example. Example 2 (Non-valid loos k) Assume that we are looking for the best slope of a linear function through data points D={(1,1), (2,5.5), (3,3)}R 2yielding the unconstrained convex parameter estimation problem min p∈[0,10] (xd,yd)∈D (p·xd−yd)2.(1) Assume further that the reduced dataset at B&B iteration kis given by Dk={(1,1), (3,3)}. Due to the convexity of (1), we can solve it also as the lower bounding problem. We observe lred k=0 at optimal solution point plb,red k=1, meaning that we can fit data points Dkexactly by the identity. We obtain loos k=(1·2−5.5)2=12.25, and lfull k=8.75 at optimal solution point plb,red k=1.5. Hence, out-of-sample lower bound loos kis larger than the full lower bound lfull k. With the help of both reduced and out-of-sample lower bound, we can enclose the full lower bound. Lemma 2 Let gcv(·; xd,yd)be Lipschitz continuous in domain Pwith Lipschitz constants Ld>0for any d =1,...,|D|.If plb,red kis feasible for (LBP), then lred k+ loos k−L·||plb,full k−plb,red k|| ≤ lfull k≤lred k+ loos k with L := |D| d=1Ld. Proof For the first inequality, we observe  (xd,yd)∈D gcv(plb,full k;xd,yd) = (xd,yd)∈Dgcv(plb,red k;xd,yd)+gcv(plb,full k;xd,yd)−gcv(plb,red k;xd,yd) = (xd,yd)∈Dk gcv(plb,red k;xd,yd)+ (xd,yd)∈D\Dk gcv(plb,red k;xd,yd) + (xd,yd)∈Dgcv(plb,full k;xd,yd)−gcv(plb,red k;xd,yd).(2) 123 622 Journal of Global Optimization (2025) 92:615–642 Since plb,red kis feasible for (LBP)and plb,full kminimizes (LBP), we have  (xd,yd)∈D gcv(plb,full k;xd,yd)≤ (xd,yd)∈D gcv(plb,red k;xd,yd) and, together with Lipschitz continuity, − (xd,yd)∈Dgcv(plb,full k;xd,yd)−gcv(plb,red k;xd,yd) =|  (xd,yd)∈Dgcv(plb,full k;xd,yd)−gcv(plb,red k;xd,yd)| additivity ≤ (xd,yd)∈D |gcv(plb,full k;xd,yd)−gcv(plb,red k;xd,yd)| Lipschitz ≤ (xd,yd)∈D Ld·||plb,full k−plb,red k||.(3) Inserting (3)in(2)gives  (xd,yd)∈D gcv(plb,full k;xd,yd)≥ (xd,yd)∈Dk gcv(plb,red k;xd,yd)+ (xd,yd)∈D\Dk gcv(plb,red k;xd,yd) −L·||plb,full k−plb,red k|| which is the same as the first inequality of Lemma 2. For the second inequality, we observe lfull k= (xd,yd)∈D gcv(plb,full k;xd,yd)≤ (xd,yd)∈D gcv(plb,red k;xd,yd) due to the optimality of plb,full kand the feasibility of plb,red kfor the full lower bounding problem (LBP). Moreover, we have  (xd,yd)∈D gcv(plb,red k;xd,yd)= (xd,yd)∈Dk gcv(plb,red k;xd,yd) + (xd,yd)∈D\Dk gcv(plb,red k;xd,yd) which concludes the proof.  Regarding the feasibility of point plb,red kfor (LBP), we note that the data-independent constraints  hcv are invariant to the dataset used. Thus, the feasibility is naturally given for any model without data-dependent constraints hcv(xd,yd;·). Based on Lemma 2,wedefinetheheuristiccombined lower bound  lcombi k:= lred k+ loos k. Note that the positivity of the cumulated Lipschitz constant Lis crucial for bounding the full lower bound below based on the combined lower bound. Due to the exhaustiveness of branching, we have diam(Pk)→0 for B&B iteration k→∞and, thus, || plb,full k− plb,red k||→0(k→∞). Even at the beginning of the B&B algorithm, i.e., for small k,the optimal points plb,full kand plb,red koften coincided in numerical experiments performed with the models investigated in Section 5. Consequently, we expect numerical advantages when 123 Journal of Global Optimization (2025) 92:615–642 629 Then, Var ⎡ ⎣gcv Dk,1(plb,red,1 k)+gcv Dk,2(plb,red,2 k) 2⎤ ⎦ = 1+Corr gcv Dk,1(plb,red,1 k), gcv Dk,2(plb,red,2 k) 2·Var gcv Dk,1(plb,red,1 k). The statement of Lemma 5follows from Dk,1and Dk,2having the same distribution as well as the scaling and additive properties of variances. As Corr [·,·]∈[−1,1]by the Cauchy-Schwarz inequality, Lemma 5implies that bagging can only improve the variance Var ⎡ ⎣gcv Dk,1(plb,red,1 k)+gcv Dk,2(plb,red,2 k) 2⎤ ⎦≤Var gcv Dk,1(plb,red,1 k). Equality can only be attained for perfectly correlated datasets giving Corr gcv Dk,1(plb,red,1 k), gcv Dk,2(plb,red,2 k)=1. Contrarily, for uncorrelated datasets, i.e., a correlation of 0, we can half the variance by picking a second subsample Var ⎡ ⎣gcv Dk,1(plb,red,1 k)+gcv Dk,2(plb,red,2 k) 2⎤ ⎦=1 2Var gcv Dk,1(plb,red,1 k). The correlation equals zero for independent random variables. However, we cannot expect different reduced datasets to be independent, since (i)there may be systematic error affecting a subset of data points similarly and (ii)the datasets are drawn from the same sample Dand may therefore contain common data points. The extent of (i)highly depends on the model and the actual dataset Dwhich are both part of the fixed problem formulation in our setup. For (ii), we have higher chances that different subsamples Dk,1and Dk,2are not intersecting, if the size of the reduced datasets Dkis small compared to the full dataset D. By default settings, we pick 10% of the data points from the full dataset for the initial dataset and add another 25% of the data points when augmenting. This means the second smallest dataset contains already a comparatively large part of the full dataset, namely 35% of all data points. For example, assume a full dataset with |D|=100. The probability that two subsamples with |Dk,∗|=10 do not intersect is about 0.33, while it is about 2.75 ×10−9for subsamples with |Dk,∗|=35. Thus, we use the proposed resampling heuristic only to update the reduced lower bound calculated based on the initial dataset. 5 Numerical results In this section, we study the computational performance of the B&B algorithm with growing datasets using the standard B&B algorithm as a benchmark. In detail, we investigate the deterministic approach from [10] as well as the heuristic approaches from Sections 3.2 and 3.3 using different augmentation rules. We run each of the approaches with and without resampling the initial dataset, cf. Section 4. Both the standard B&B algorithm and all discussed algorithmic approaches of the B&B algorithm with growing datasets are available in our open-source solver MAiNGO v0.8.2. 123 630 Journal of Global Optimization (2025) 92:615–642 Table 2 Overview on general properties and references of the models and data used in the case study Name n|D|dim(x)Data Opt. class Original references Model Data Process systems engineering EOS262 9 262 2 synthetic, exact MINLP [21,22][10] EOS2262 9 2262 2 synthetic, exact MINLP — ”—— ”— EOS262noisy 9 262 2 synthetic, noisy MINLP — ”—— ”— EOS2262noisy 9 2262 2 synthetic, noisy MINLP — ”—— ”— IHMcon 20 284 1 measured NLP [23][24,25] IHMunc 10 284 1 measured NLP — ”—— ”— kinetics 5 446 1 measured DO [1,2,26][2] EIS 7 26 1ameasured NLP [27][28,29] Biochemistry TSP 12 20 8 synthetic, exact DO [30,31][10] TSPnoisy 12 20 8 synthetic, noisy DO — ”—— ”— Machine learning GMMcon 6 272 1 measured NLP [32][33,34] GMMineq 5 272 1 measured NLP — ”—— ”— trainANN 21 220 3 measured NLP [35,36][37,38] aThe model output is a complex number, i.e., dim(y)=2 Benchmark libraries from mathematical programming, e.g., MINLPlib [39], PrincetonLib [40], and the COCONUT benchmark [41], only account for parameter estimation problems (PE) considering small datasets which do not require our extension. Besides, our focus is the solution of real-world applications extending the numerical proof-of-concept obtained in [10]. Hence, we collected different models from both literature and our previous work stemming from process systems engineering, machine learning, and biochemistry. These models cover different classes of optimization problems including mixed-integer nonlinear programs (MINLPs), nonlinear programs (NLPs), and dynamic optimization problems (DOs) with up to n=21 unknown parameters and a different number and dimension of data points in full dataset D,seeTable2. The exact mathematical expressions are provided in Online Resource 1. For ready-to-use implementations of the model as well as the exact data points used refer to the aforementioned repository GloPSE2. All computations are performed on Intel Xeon Platinum 8160 processors “SkyLake” (frequency 2.1GHz, RAM =3.75GB for a single node). Even small deviations in computational times for processing one node may significantly change the number of nodes processed within our CPU time limit of 23h when adding up. Consequentially, the final lower and upper bound may change as well. Thus, we repeat each run 5 times and report the median of the final bounds and, in case of convergence, CPU times. Note that we choose a CPU time limit of 23h such that we can guarantee a total runtime of 24h including the initialization of the model and preparations of the output data. We set the relative optimality tolerance to εR=0.1. The absolute tolerance is fixed to εA=0.01 for all approaches using (SSE). When using (MSE), the objective is scaled by the number of data points. In analogy, we scale the absolute optimality tolerance to model specific values εA=0.01/|D|when using (MSE). In each optimization run, we use Ipopt version 3.12.12 [42] for running 3 local searches as a pre-processing step in the root node. Moreover, 123 Journal of Global Optimization (2025) 92:615–642 631 Fig. 5 Lower bounds of nodes in Nchanged before (◦) and after (×) post-processing and change of final lower bound ldue to post-processing compared to final upper bound uand relative optimality tolerance εR we use McCormick relaxations [43,44] calculated with MC++ [45] and a linearization in the midpoint of the parameter intervals to obtain a linear program (LP) for lower bounding. For models IHMcon, IHMunc, EIS, TSP, TSPnoisy, GMMcon, GMMineq, and trainANN, we solve the LP with linear optimizer CLP v1.17.0 [46]. As the resulting LP of the EOS models with noisy measurement data as well as model kinetics seems to be numerically difficult for CLP and alternative LP solver CPLEX [47], we use interval extensions [48] calculated with FILIB++ [49] for models EOS262, EOS2262, EOS262noisy, EOS2262noisy, and kinetics. When using pure interval extensions, we disable the optimality-based bound tightening [50]. The upper bounding problem is solved with local optimizer LBFGS [51,52] implemented in the NLOPT toolbox v2.5.0 [53] for all models. To minimize deviations caused by small differences in the CPU time required for processing a node, we limit the number of steps performed by local solvers Ipopt and LBFGS rather than the CPU time used for each local search. The complete listing of settings is provided in Online Resource 2. 5.1 Evaluation of post-processing At first, we study the post-processing step to check whether we can expect the heuristic approaches to converge to the global solutions. Table 3summarizes the number of nodes tracked for post-processing, the implications for the lower bounds reported, and the maximum CPU time used for post-processing. Note that these are the statistics of run 1 out of the 5 repetitive runs only. The number of nodes tracked as well as which specific nodes are tracked may vary if the number of processed nodes differs. However, we expect to make the largest mistakes early, namely with small datasets, which is covered by all 5 repetitive runs. More importantly, note that the statistics of models IHMcon, IHMunc, EIS, and trainANN are not listed, just as the statistics of the MSE heuristic for models TSP and TSPnoisy, since no nodes are tracked for post-processing in these cases. In these cases, we do not prune based on a reduced dataset. This means we have a deterministic procedure so far at the cost of a large B&B tree. Similarly, in runs with |Npostpro|<100 it seems to be hard to prune based on the reduced datasets. As an exception, |Npostpro|<100 comprises also the case where we need to process less than 100 nodes for convergence, cf. TSP model in Table 3. We observe that for only 4 out of 13 models there are nodes where the lower bound lfull k calculated in post-processing falls below the final upper bound, i.e., where we may have made a wrong pruning decision due to data reduction. Out of these, the final lower bound is updated solely for 3 models, namely when using the SSE heuristic applying augmentation rules scaling and tol for model EOS262 and augmentation rule tol for models GMMcon and GMMineq. We pick two exemplary cases to depict the changes due to post-processing in Figure 5.In both cases, the nodes which affect the final lower bound  lhave been only just pruned with  lcombi k(1−εR)·u. In model EOS262, the changes due to post-processing are essentially nonexistent: we find numerically insignificant differences of 10−14 to 10−18. In contrast to that, the final lower bound of models GMMcon and GMMineq change significantly, e.g., 123 632 Journal of Global Optimization (2025) 92:615–642 Table 3 Post-processing statistics for the models with tracked nodes, i.e., with Npostpro =∅, including the number of nodes tracked, the number of nodes where the lower bound changed significantly within post-processing, a flag indicating whether the final lower bound was changed as well, and the maximum CPU time over all augmentation rules Model Re- |Npostpro|No. of nodes wrongly prunedaFinal lower bound changed?bMaximum sampling? CPU time [s] SSE heuristic using augmentation rules CONST /SCALING /TOL EOS262 No 100/ 100 / 100 0 / 2/2 0/ 1 / 1 0.7 Yes 100/ 100 /100 0/ 3 / 3 0/ 1/1 0.9 EOS2262 no 100/ 100 / 100 0 / 0/0 0/ 0 / 0 8.2 Yes 100/ 100 /100 0/ 0 / 0 0/ 0/0 9.9 EOS262noisy no 22/ 72 / 72 0/ 0 / 0 0/ 0 /0 0.6 Yes 22/72/72 0/0/0 0/0/0 1.0 EOS2262noisy no 52/ 100 / 100 0/ 0 / 0 0/ 0/0 7.1 Yes 66/ 100 /100 0/ 0 / 0 0/ 0/0 7.9 kinetics no 100/ 100 / 100 0 / 0/0 0/ 0 / 0 1.6 Yes 100/ 100 /100 0/ 0 / 0 0/ 0/0 1.6 TSP no 100/ 100 / 4 0 / 0/0 0/ 0 / 0 1.3 Yes 100/ 100 /4 0/ 0 / 0 0/ 0/0 0.3 TSPnoisy no 16/ 100 / 0 0 / 0/0 0/ 0 / 0 0.2 Yes 0/100/0 0/0/0 0/0/0 0.2 GMMcon no 44/ 100 / 100 0 / 6/74 0/ 0 / 1 0.1 yes 44/ 100 / 100 0 / 2/70 0/ 0 / 1 0.1 123 Journal of Global Optimization (2025) 92:615–642 633 Table 3 continued Model Re- |Npostpro|No. of nodes wrongly prunedaFinal lower bound changed?bMaximum sampling? CPU time [s] GMMineq no 100/ 100 / 100 0 / 0/96 0/ 0 / 1 0.2 yes 100/ 100 / 100 0 / 0/99 0/ 0 / 1 0.2 MSE heuristic using augmentation rules CONST /OOS /TOL EOS262 no 100/ 100 / 100 0 / 1/0 0/ 0 / 0 0.5 yes 100/ 100 / 100 0 / 0/0 0/ 0 / 0 1.0 EOS2262 no 100/ 100 / 100 3 / 1/1 0/ 0 / 0 6.4 yes 100/ 100 / 100 2 / 4/1 0/ 0 / 0 7.1 EOS262noisy no 12/100 / 100 0 / 0 /0 0/0 / 0 0.6 yes 30/12/92 0/0/0 0/0/0 0.8 EOS2262noisy no 4/ 76 / 76 0/ 0 / 0 0/ 0 /0 5.2 yes 12/76/76 0/0/0 0/0/0 5.8 kinetics no 100/ 100 / 100 0 / 0/0 0/ 0 / 0 1.4 yes 100/ 100 / 100 0 / 0/0 0/ 0 / 0 1.7 GMMcon no 100/ 100 / 100 9 / 17/17 0/ 0 / 0 0.1 yes 100/ 76 / 100 9 / 12/8 0/ 0 / 0 0.1 GMMineq no 100/ 100 / 100 1 / 4/4 0/ 0 / 0 0.1 yes 100/ 100 / 100 1 / 4/4 0/ 0 / 0 0.0 aNodes with lfull k<u, see Section 3.2 b0 = no, 1 = yes 123 634 Journal of Global Optimization (2025) 92:615–642 from 723.7 to 406.7 for GMMineq when using resampling, extending the relative optimality gap from 15% to 86%. For 95% of 156 cases, including the 2 heuristic approaches with 3 different augmentation rules each for 13 models with and without resampling, the final lower bound was not changed and only in 4 cases the change was significant. In particular, the final lower bound was changed only when using the SSE heuristic. Note that the post-processing is very cheap: over all 5 repetitions of the case study we measured a maximum runtime of 11.3s. In conclusion, even the heuristic approaches converge towards the global solution in most of the cases. The post-processing procedure detects the remaining cases with small computational effort and provides means for deciding whether the heuristic pruning significantly distorts the solution obtained. 5.2 Comparison of computational performance Finally, we study the performance of the B&B algorithm with growing datasets with the standard B&B algorithm as a benchmark. Only 4 models converge within the CPU time limit of 23h for at least one of the approaches including all 3 models with exact data, namely EOS262, EOS2262, and TSP, compare Table 4. For models EOS262 and EOS2262, the deterministic approach is the fastest with decreasing the runtime of the standard B&B algorithm by a factor of 3 and up to 4, respectively. For the TSP model, the MSE heuristics finds the global solution in the root node, resulting in a runtime of 2s, while the standard B&B algorithm takes about 2h for convergence with (SSE) and hits the CPU time limit with (MSE). Note that for exact data, the final lower bound equals the natural lower bound of 0 except for numerical tolerances. Thus, using reduced datasets allows for reducing the computational costs while retaining the tightness of the bounds on the optimal solution. Contrarily, the standard B&B algorithm performs best for the EIS model: it converges within 8h, whereas all approaches of the B&B algorithm with growing datasets hit the CPU time limit. In the EIS model, we fit a model with 7 unknown parameters to a full dataset D R×R2with |D|=26 data points and, with the default settings, to reduced datasets containing 3, 10, 17, and 24 data points. Consequently, this parameter estimation problem is prone to overfitting, in particular, when using reduced datasets. In other words, the computational savings due to the data reduction come at the cost of looser lower bounds. Out-of-sample estimation allows for improving upon the quality of the lower bounds calculated based on reduced datasets, cf. Figure 6. Figure 6provides an overview of the range of final lower and upper bounds obtained by the different algorithms within the CPU time limit. Note that the lower bounds include, where applicable, the correction performed in the post-processing step and that the exact values are given in Online Resource 1. The first row of Figure 6contains the models considering exact data such that the final lower bound is equal to zero except for numerical tolerances and the optimal solution is in the range of the optimality tolerance. For the remaining models, all approaches give a similar or even the same solution after the CPU time limit, compare the upper bounds depicted in Figure 6. For 7 out of 10 models considering noisy data, the SSE heuristic gives the best, i.e., largest, final lower bound. For 2 of these 10 models, the SSE heuristic and the standard B&B algorithm perform similarly. Only for IHMcon, the standard B&B algorithm finds the best lower bound. In accordance with the results of Section 5.1, we therefore conclude that the SSE heuristics allows to significantly increase the final lower bound within a given CPU time limit. 123 Journal of Global Optimization (2025) 92:615–642 635 Table 4 Total CPU time needed for convergence of the B&B algorithm. Only models which converge within the CPU time limit for at least one of the algorithmic approaches, namely the standard B&B algorithm with objective (SSE)and(MSE), or the deterministic and heuristic approaches of the B&B algorithm with growing datasets using different augmentation rules are listed. Bold numbers indicate the fastest runtime for the respective model Standard B&B Deterministic Appr. (SSE) SSE heuristic MSE heuristic Model (SSE)(MSE)const scaling combi const scaling tol const oos tol Without resampling EOS262 2.0h 2.1h 1.7h 44min 40 min 13.3h –a–a7.9h –a5.8h EOS2262 –a–a–a7.5h7.5h– a–a–a–a–a–a EIS 9.2h 8.0h– a–a–a–a–a–a–a–a–a TSP 1.9h –a–a40min –a–a–a9.2h 2s2s2s With resampling EOS262 2.0h 2.1h 1.8h 47min 38 min 12.8h –a–a12.1h –a5.7h EOS2262 –a–a–a5.5h6.0h –a–a–a–a–a–a EIS 9.2h 8.0h– a–a–a–a–a–a–a–a–a TSP 1.9h –a–a1.0h –a–a–a7.1h 2s2s2s aHitting CPU time limit of 23h 123 636 Journal of Global Optimization (2025) 92:615–642 Fig. 6 Final lower bounds lkor  lk(strictly below dashed line) and final upper bounds u(above dashed line) on the optimal solution of the models for any of the approaches, namely the standard B&B algorithm (std) with objectives (SSE)and(MSE) as well as the deterministic approach (det), the SSE heuristic (heuS) and the MSE heuristic (heuM) of the B&B algorithm with growing datasets using the augmentation rules depicted in Figures 3and 4. If applicable, the final lower bounds are subject to post-processing. The squares in the logarithmic plots indicate lower bounds equal to 0 We recall that models IHMunc and IHMcon are mathematically equivalent. While the optimization variables are coupled via equality constraints in the former, we use this relation as an explicit evaluation function for the latter, compare also the full-space and reducedspace formulations discussed in [2,54]. While the final lower bound reported for IHMunc differs from 0 only when using the heuristic approaches, only the standard B&B algorithm and the deterministic approach determine lower bounds larger than 0 for IHMcon. The best lower bound for model IHMunc is in order of magnitude 106, while the best lower bound for model IHMcon is in order of magnitude 105. Thus, applying the heuristic approaches of the 123 Journal of Global Optimization (2025) 92:615–642 637 B&B algorithm with growing datasets to the reduced-space formulation performs best for the IHM model. In analogy, the equality constraint in GMMcon is used to fix an unknown parameter in GMMineq. In this case, the reduced-space formulation GMMineq seems to be slightly advantageous for finding a good upper bound, while the lower bounds obtained for models GMMcon and GMMineq are similar for the different approaches. When comparing the CPU times with and without resampling, cf. Table 4, there is no general tendency. On the one hand, resampling requires more computational resources. On the other hand, resampling affects the lower bound. Since MAiNGO uses the solution point of the lower bounding problem to initialize the local solver for the upper bounding problem, we may obtain a much better upper bound with resampling by chance. In this case, we need significantly less iterations for convergence. Apart from this, resampling does not seem to affect the final lower bounds in most of the cases, see Figure 6. In fact, the deviation in the number of nodes processed over the 5 repetitive runs affects the results in some cases as much as the choice whether to use resampling, e.g., the final lower bounds reported for model EOS2262. Noting that the resampling heuristic was introduced to reduce the variances of the reduced lower bounds, we infer that the B&B algorithm with growing datasets does not suffer from large variances. 6 Conclusions and outlook We investigate out-of-sample estimation for enhancing the B&B algorithm with growing datasets proposed in our previous work [10]. In detail, we combine the lower bound calculated in a B&B node based on a reduced dataset with an out-of-sample evaluation obtaining the combined lower bound. Although the combined lower bound is only a heuristic lower bound, we expect it to be close to the lower bound calculated based on the full dataset. We use this heuristic lower bound for extending the deterministic approach presented in [10]aswellas for pruning which results in two heuristic approaches, the so-called SSE and MSE heuristic. To detect and quantify potential mistakes made by the SSE and MSE heuristic, we introduce a post-processing check of the final lower bound after the B&B algorithm terminated. We compared the performance of the different approaches of the B&B algorithm with growing datasets and the standard B&B algorithm based on 13 real-world applications from both literature and our previous works. For the estimation problems with exact data, the deterministic approach of the B&B algorithm with growing datasets yields the fastest runtimes. Most of our problems with error-prone data do not converge within the given CPU time limit. For these models, all approaches, including the standard B&B algorithm, find similar upper bounds. The SSE heuristics yields the best, i.e., largest, final lower bound for 70% of the models with noisy data and a similar value as the standard B&B algorithm for 20% of these models. In turn, the SSE heuristic may introduce an error into the pruning procedure. However, our results suggest that both the SSE and the MSE heuristic are almost deterministically converging towards the global solution, where the proposed post-processing procedure allows to detect and quantify the exceptions from this assumption. Apart from this, we show in theory that we can likely decrease the variance of the lower bounds calculated based on reduced datasets when resampling the dataset. However, the actual numerical performance of the B&B algorithm with growing datasets for the models of our case study is hardly affected by the use of the resampling heuristic. Since we always fit a comparatively small number of unknown parameters to a large dataset, we seem to obtain small variances for the reduced lower bounds making resampling techniques like bagging less 123 638 Journal of Global Optimization (2025) 92:615–642 advantageous, cf. [20]and[15, Section 8.2]. Note that the good generalization performance of our solution of model trainANN, see Online Resource 1, is another indication for having small variances. The computational advantage of the B&B algorithm with growing datasets depends on the choice of the reduced datasets. On the one hand, estimation results based on a reduced dataset may be distorted by overfitting or even identifiability issues if the full dataset is already small compared to the number of unknown parameters. We have shown that applying outof-sample estimation allows for tightening the lower bound calculated based on the reduced dataset. However, if the data reduction is insignificant in absolute numbers, the resulting computational savings cannot compensate for the remaining gap to the full lower bound. As a consequence, the B&B algorithm with growing datasets is still outperformed by the standard algorithm in such cases with small datasets. On the other hand, we expect that the performance of the B&B algorithm with growing datasets profits from exploiting knowledge about the dataset for the data reduction. For example, the reduced lower bound may be a better fit to the full lower bound if the reduced datasets contain a sufficient amount of measurements with the largest measurement errors. An extension allowing for user-given reduced datasets merits therefore careful attention. We make the parameter estimation problems of our case study openly accessible via repository GloPSE2. In future work, we aim at extending the case study with estimation problems containing data-dependent constraints, e.g., fitting binary fluid systems with constraints on the measured mole fractions [55] and the optimization of Gaussian processes with constraints on each of the data points of the training set [56]. For this, we aim at implementing a specific treatment of data-dependent constraints to allow for an efficient handling of reduced datasets for these models. For all computations, we used our open-source solver MAiNGO which uses McCormick relaxations [43,44] to calculate convex underestimators for the lower bounding problem. MAiNGO can therefore efficiently handle reduced-space formulations [2,54], meaning that the size of the dataset does not affect the number of optimization variables but solely the computational effort for function evaluations. In contrast to that, many DGO solvers like BARON [57,58], ANTIGONE [59], and SCIP [60,61] use the auxiliary values method (AVM) [6, 62–64] for obtaining convex underestimators. An implementation of the B&B algorihm with growing datasets in these AVM-based solvers is of high interest since the AVM method may add an auxiliary optimization variable for each of the data points and we therefore expect AVM-based solvers to profit even stronger from the use of growing datasets. Supplementary Information The online version contains supplementary material available at https://doi. org/10.1007/s10898-025-01514-4. Acknowledgements This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant MI 1851/10-1 “Parameter estimation with (almost) deterministic global optimization”. Simulations were performed with computing resources granted by RWTH Aachen University under project rwth1563. Susanne Sass is grateful for her association to the International Research Training Group (DFG) IRTG-2379 “Hierarchical and Hybrid Approaches in Modern Inverse Problems” under grant 333849990/GRK2379. We thank Dominik Bongartz for his support in the implementation of the novel approaches within MAiNGO. Moreover, we thank Ian H. Bell, Tanvir Rahat, Alexander Echtermeyer, and J. Raphael Seidenberg for their assistance with setting up the EOS, TSP, IHM, and EIS model, respectively. We are grateful to Ian H. Bell for providing the data for the EOS model as well as to Niklas Thissen, Stefanie Khan, and Anna K. Mechler for providing experimental measurements for the EIS model. Moreover, we thank the anonymous reviewer whose comments helped us with improving and clarifying this article. Author Contributions Conceptualization: Susanne Sass, Angelos Tsoukalas; Methodology: Susanne Sass, Alexander Mitsos, Nikolay I. Nikolov (particularly of Section 4), Angelos Tsoukalas; Software: Susanne Sass; 123