Solving linear multiplicative programs via branch-and-bound: a computational experience

In this paper, linear multiplicative programs are approached with a branch-and-bound scheme and a detailed computational study is provided. Several underestimation functions are analyzed and various partitioning criteria are presented. A particular class of linear multiplicative programs, useful to solve some applicative bilevel problems, is considered from a theoretical point of view to emphasize an efficient solution method. Detailed results of the computational study are provided to point out the performances provided by using various underestimation functions and partitioning criteria, thus improving some of the results of the current literature.


Introduction
In the literature, linear multiplicative programs have been widely studied due to their importance from both theoretical and applicative point of views.These problems, strictly related to quadratic programming and bilinear programming, are used in plant layout design, portfolio and financial optimization, VLSI chip design, robust optimization, network flows (see, e.g., Cambini and Martein 2009;Cambini and Sodini 2008;Cambini and Salvi 2010;Gupta 1995;Horst and Pardalos 1995;Horst 38 Page 2 of 32 and Tuy 1996;Horst et al. 2001;Konno and Kuno 1992;McCarl et al. 1977;Mjelde 1983;Ryoo and Sahinidis 2003;Tuy 2016 and references therein).
From the computational point of view, various approaches have been proposed.For instance, Wang et al. (2012) proposed a branch-and-bound algorithm for global minimization of a generalized linear multiplicative programming; Jiao et al. (2012) proposed a branch and bound algorithm based on the computation of subsequent solutions of the series of linear relaxation programming problems; in a very recent paper, Jiao et al. (2023) proposed a branch-reduction-bound algorithm, based on outer space search and branch-and-bound framework.
A second field of literature develops methods based on an eigenvectors approach.The eigenvectors approach has been widely used despite actually it has some drawbacks.It is well known that the eigen-decomposition of a quadratic function is an heavy task from both a computational and a numerical point of view and it does not care about the particular structure of linear multiplicative functions.For all these reasons, in the recent years, various papers have introduced new procedures developed without the use of eigen-decompositions (see, e.g., Shen et al. 2020Shen et al. , 2022;;Wang et al. 2012;Zhou et al. 2015).
In this paper, various underestimation functions to be used in the branch and bound procedure are introduced and discussed both with the eigenvectors approach and without it.Then, a full computational test is provided to highlight the best performing functions.
In addition, a motivating example is presented in the last part of the paper.Special structures of multiplicative problems, in facts, arise in several bilevel programming problems of leader-follower type (see Dempe 2020, for example).Some underestimation functions introduced in the paper can be used to improve the algorithmic procedure adopted to tackle this class of problems.
The main contributions of this paper are: • To describe a unified framework for dealing with linear multiplicative programs from both a theoretical and computational point of view; • To propose various underestimation functions to be used in the branch and bound procedure: quadratic, linear, difference of two convex functions (D.C. functions), and eigenvectors based underestimation functions; • To perform a detailed computational test that compares the different underestimation functions under various partition methods to identify the most promising ones; • To characterize a special case of multiplicative programs strictly related to bilevel optimization and to define the more appropriate underestimation function to solve them.
The paper is organized as follows.In Sect.2, the main definitions and preliminary results are given.In addition, the criteria for the splitting process of the branchand-bound approach are analyzed.On this basis, Sect. 3 is devoted to the study of quadratic, eigenvectors based and linear underestimation functions.Then, in Sect.4, the detailed results of a wide computational experience are provided and fully discussed, giving a detailed view of the computational aspects of the solution method and improving some of the results of the current literature.Furthermore, Sect. 5 points out the behavior of the proposed underestimation functions in a particular class of linear multiplicative programs very useful in applicative bilevel programming.Finally, a section with the conclusions is given.

Definitions and preliminary results
The aim of this section is to define the problem and provide the main preliminary results which will allow the development of the paper.In this light, firstly the problem is defined and then the concept and properties of underestimation functions are given.On this basis, a detailed description of a branch-and-bound scheme to solve the problem is provided, in order to approach it in a unifying framework with respect to different underestimation functions and different splitting criteria.Finally, the quadratic form associated to a linear multiplicative function is recalled as well as the eigenvector-based decomposition of such a quadratic form.

Definition of the problem
From now on, ℝ will denote the set of real numbers while ℝ = ℝ ∪ {−∞, +∞} will be the affinely extended real number system.
Definition 1 Let P be the following minimization problem: where S ⊆ ℝ n is a nonempty polyhedron defined as with A in ∈ ℝ m×n , b in ∈ ℝ m , A eq ∈ ℝ r×n , b eq ∈ ℝ r , and l b , u b ∈ ℝ n , while f ∶ ℝ n → ℝ is a linear multiplicative function defined as

Underestimation functions
The use of suitable underestimation functions of f is needed to solve Problem P by means of a branch-and-bound approach.moreover, Er ∶ S → ℝ is the corresponding error function defined as The following useful properties hold.
Proof For any x ∈ S , it results:

◻
The following further result will be used in the next subsection as an algorithmic stopping criterium.
Lemma 2 Let S ⊆ ℝ n be nonempty.Let Φ ∶ S → ℝ be an underestimation function of f ∶ S → ℝ , with Er ∶ S → ℝ its error function.If x ∈ arg min S Φ and Er(x) = 0 , then x ∈ arg min S f .Proof For all x ∈ S , it results: hence, the thesis follows.◻ Finally, it is worth noticing that some underestimation functions of f will be obtained by first rewriting f in D.C. form.Let us recall that a function f ∶ ℝ n → ℝ is said to be in D.C. form if it is expressed as f (x) = q 1 (x) − q 2 (x) , with q 1 ∶ ℝ n → ℝ and q 2 ∶ ℝ n → ℝ convex functions (see, e.g., Cambini and Salvi 2009, 2010 and  references therein).
In this light, the following result will be useful to deduce underestimation functions of f in the case f is rewritten in D.C. quadratic form.
Then, for all y ∈ [y L , y U ] , the error given by the use of −(y Proof For all y ∈ [y L , y U ] , the error given by the use of −(y and the thesis follows.◻

A branch-and-bound scheme
In order to solve Problem P by using a branch-and-bound approach (see, e.g., Bajaj and Faruque Hasan 2020;Cambini andSodini 2005, 2008;Cambini andSalvi 2009, 2010;Fampa et al. 2017;Gerard et al. 2017;Shen et al. 2020 and references therein), the following operative scheme will be considered: • the feasible region will be iteratively partitioned in multidimensional rectangles, • the function f will be "relaxed" over the single partitions by means of a convex underestimation function Φ, • the convex "relaxed" subproblems will be solved • the feasible solution having the smallest value of f will be maintained.
Specifically speaking, such a branch-and-bound approach will be implemented by means of: • a Priority Queue (PQ) used to store the single partitions sorted with respect to the Lower Bound (LB) found minimizing the convex "relaxed" subproblems over the partitions; • a feasible point Sol and a value UB corresponding, respectively, to the incumbent best feasible solution found and its image UB = f (Sol); • a convex underestimation function Φ needed to solve the subproblems over the vari- ous partitions.
The priority queue PQ is used to speed up the solution method (at the cost of memory usage of course) since the partition with smaller LB is always known and since it does not need of any periodic "pruning" process (in the "pruning" process of a branch-andbound scheme, the stored partitions having a LB not smaller than UB are cancelled since they cannot improve the incumbent best feasible solution).
The following commands are aimed to describe the way the priority queue PQ is managed: The command PQisempty tells whether the PQ is empty or not; the command PQsmallest provides the smallest LB of the partitions stored in the PQ; the command PQadd adds to the PQ a partition as well as the corresponding Optimal Solution (opt) and its optimal value LB = Φ(opt) obtained minimizing Φ over the partition itself.Clearly, PQ is such that smaller is the value of the LB, the higher the priority of the partition in the PQ will be.The last command PQextract removes from the PQ the partition having the smallest LB and provides as the output all the data stored.
Remark 1 If the smallest LB is such that LB ≥ UB , then the partitions in the PQ are not able to improve the incumbent best feasible solution Sol and hence PQ can be emptied (being PQ a priority queue, the "pruning" process becomes just a final stopping condition and it is no more periodic).
The following subprocedure "CheckPartition()" minimizes the underestimation function Φ over a given partition Π , improves the value of UB and stores the data in the PQ when the potential relative improvement is such that with and relTol > 0 be a chosen tolerance.According to Remark 1, if LB ≥ UB then the partition is not added to the queue and discarded ("pruning").
The overall solution branch-and-bound scheme is described by the procedure "Solve()".
( Firstly, variables PQ, UB and Sol are initialized and the starting smallest partition Π 0 containing S is found and checked.Then, the iterative phase starts and continues up to either the PQ is emptied or the stored partitions cannot improve anymore UB ("pruning").At each iteration, the partition with the smaller LB is extracted from the PQ, a branching variable is selected, and the partition split accordingly (multiway branching has been shown to provide poor results, see for example Gerard et al. 2017).The two new partitions are then checked.Finally, at the end of the iterative process, outputs are set.

Remark 2
The value of UB is fundamental to improving the performance of the algorithm since it is used to discard the not useful partitions.For this reason, the feasible points found while looking for Π 0 should be used to improve values UB and Sol.
Notice that the convergence of the proposed method has been widely discussed in the literature (see, e.g., Bajaj and Faruque Hasan 2020;Cambini andSalvi 2009, 2010;Fampa et al. 2017;Gerard et al. 2017;Shen et al. 2020 and references therein).Specifically speaking, since the partitions will be split with respect to values not "close" to its boundaries (see Sect. 4), then the tolerance parameter relTol > 0 guarantees that condition (1) in subprocedure "CheckPar- tition()" will become false after a sufficiently large number of iterations.The correctness of the method follows since just feasible solutions are evaluated to improve the incumbent best solution and since the whole feasible region is analyzed.This is known to be an NP-Hard problem and in the worst case many local but not global optimal solutions can be found.Some further choices are finally needed to complete the description of the solution process: • Which underestimation function Φ(x) should be used?tight underestimation functions improve the algorithm performance, moreover the underestimation function determines the set of branching variables; • Which branching variable should be chosen to split the current partition?
• With respect to which value of the branching variable the current partition should be split?
Actually, another fundamental choice has been already made: • At each iteration, the partition with the smaller LB is selected and analyzed.
This criterium is aimed to look for feasible solutions having small values, thus allowing to improve UB as much as possible and to increase as much as possible the number of partitions discarded by means of the stopping "pruning" condition.

A raw approach
Problem P is a particular quadratic (usually indefinite) program since f(x) can be rewritten as: with Quadratic indefinite programs can be efficiently solved with a branch-and-bound approach by means of a suitable eigenvectors-based decomposition of the objective function (see, e.g., Cambini andSodini 2005, 2008;Fampa et al. 2017).Specifically speaking, since Q is a symmetric matrix, there exists an orthonormal matrix The diagonal elements of D are the eigenvalues 1 , … , n ∈ ℝ of Q , while the orthonormal columns u 1 , … , u n ∈ ℝ n of U are the corresponding eigenvectors of Q .
As a consequence, it results: Thus, by means of the sets of indices and the vectors the quadratic component of f can be rewritten as follows: In this way, f can be expressed in the following D.C. form: The following branching variables are suggested by (3) for all i ∈ Λ − : so that Moreover, for all With respect to branching variables , rectan- gle , is the smallest partition Π 0 containing S. The following underestimation function can then be stated by means of (3).
Theorem 1 Let f ∶ ℝ n → ℝ be expressed as in (3).Then, the following convex quadratic function is an underestimation function for f over S ∩ L , U : and for all i ∈ Λ − , from Lemma 3 it yields: Hence, Φ 0 (x) follows trivially from (3).Moreover, it results:

Splitting process
Assume to be in the iterative process and that the partition having the smallest lower bound LB have been extracted from PQ by means of the command " [Π, opt] ∶= PQextract() ".It is now necessary to choose a branching variable in order to split the partition Π .Assume, for example, that underestimation (u 0 ) is used (results are analogous for all underestimations functions which will be proposed in the next Section), so that Π = L , U and: with maximum error value Two criteria are generally used to determine the branching variable: • The largest interval L i , U i : this criterion is aimed to reduce as much as possible the error maximum value; on the other hand, it does not use the feasible point opt = arg min x∈S∩Π Φ 0 (x) which may suggest where to look for the opti- mal solution of the problem.
• The largest error addend : this criterion is aimed to reduce as much as possible the error value at the point opt, hence tightening the underestimation as much as possible close to opt.
Once the branching variable is chosen, the question is how to split the corresponding interval L i , U i , that is, how to determine the value * i so that In this light, possible choices are: • The medium point of the interval L i , U i : the use of is aimed to reduce as much as possible the error maximum value, again the feasible point opt is not used; • The value opt i = v T i opt : in this case the optimal solution of the underestima- tion is used, but this value may be close to the boundaries L i and U i thus highly increasing the number of iterations needed to solve the problem and hence increasing the convergence time itself; • A linear combination of M i and i could be used to take into account of both the error maximum value and the solution opt.
These criteria can be summarized and linked as follows: • Blindly reduce as much as possible the error maximum value: choose the largest interval L i , U i and split it in the middle with M i (see, e.g., Cambini and Sodini 2005;Shen et al. 2020); • Use the "infos" given by opt: choose the largest error addend Cambini andSalvi 2009, 2010;Fampa et al. 2017).
The recent literature shows that the latter opportunity is the most performing one, and is the one that will be used in the computational test described in Sect. 4. Finally, notice that, by means of Lemma 2, the splitting process should be performed only if the largest error addend is grater than a suitable tolerance errTol > 0.

Specific underestimation functions
The eigenvectors-based approach described in Sect.2.4 actually has some drawbacks.First of all, the eigen-decomposition of a quadratic function is an heavy task from both a computational and a numerical point of view.Moreover, such an eigen-decomposition does not take into account of the particular structure of linear multiplicative functions.In this very light, in the recent literature various papers aimed to approach linear multiplicative problems without the use of eigendecompositions (see, e.g., Shen et al. 2020Shen et al. , 2022;;Wang et al. 2012;Zhou et al. 2015).
The aim of this section is to state some underestimation functions for f not using the eigenvectors of matrix Q , in order to efficiently solve Problem P with- out the computational and numerical troubles of eigenvectors computing.

Linear underestimation functions
In the recent literature (see Shen et al. 2020, for example) some linear underestimation functions have been used to solve linear multiplicative problems by means of a branch-and-bound approach.Usually these underestimations are not tight, for this very reason in this subsection some further linear underestimations will be studied.Recalling that: the following 2p branching variables can be considered: , and ∶= ( i ) i=1,…,p .With respect to branching variables and , rectangle , ∩ , is the smallest parti- tion Π 0 containing S.
The following linear underestimation functions can then be stated.
Theorem 2 Let f ∶ ℝ n → ℝ be expressed as in (2).Then, the following linear func- tions are underestimation functions for f over S ∩ L , U ∩ L , U : Hence, the thesis follows by means of simple calculations.
38 Page 14 of 32 1976).Notice also that in the case L ≥ 0 and L ≥ 0 , i = 1, … , p , linear underesti- mations of the kind have been used, and that these are far less tight than the one proposed in this subsection.

Quadratic underestimation functions
The aim of this subsection is to state underestimation functions of f, tighter than the linear ones, by properly rewriting f in D.C. form.
Theorem 3 Let f ∶ ℝ n → ℝ be defined as in (2).Then, f can be rewritten in the fol- lowing D.C. forms: Proof Firstly, one get: Then, by opportunely replacing each of them in (2), the thesis follows.◻ The following underestimation function can be stated by means of (i) of Theorem 3 and the 2p branching variables described in the previous subsection.
Theorem 4 Let f ∶ ℝ n → ℝ be expressed as in (i) of Theorem 3.Then, the fol- lowing convex quadratic function is an underestimation function for f over S ∩ L , U ∩ L , U : with and for all i = 1, … , p , from Lemma 3 it yields: Hence, Φ 3 (x) follows trivially from (i) of Theorem 3.Moreover, it results:

◻
In similar way, for all i = 1, … , p , the following branching variables are sug- gested by (ii) and (iii) of Theorem 3: so that Moreover, for all L ∶= ( L i ) i=1,…,p , L ∶= ( L i ) i=1,…,p ∈ ℝ p such that L ≤ U , the following rectangle is introduced: 38 Page 16 of 32 Let ∶= ( i ) i=1,…,p and ∶= ( i ) i=1,…,p .With respect to branching variables , rectangle , is the smallest partition Π 0 containing S. In this light, the following underestimation functions can be stated by means of (ii) and (iii) of Theorem 3.
Theorem 5 Let f ∶ ℝ n → ℝ be expressed as in (ii) and i(ii) of Theorem 3.Then, the following convex quadratic functions are underestimation functions for f over S ∩ L , U , respectively: and with 1 16 2 the maximum error for Er 4 (x) and Er 5 (x) , respectively, obtained at (c and for all i = 1, … , p , from Lemma 3 it yields: Hence, Φ 4 (x) and Φ 5 (x) follow, respectively, from (ii) and (iii) of Theorem 3. Moreo- ver, it results: and ◻ Remark 4 Notice that (u 4 ) dominates (u 5 ) since Er 5 (x) = 2 ⋅ Er 4 (x) .For this very rea- son, (u 5 ) has been given just for the sake of completeness and will no more be used in the rest of the paper.

Further hybrid underestimation functions
For the sake of completeness, some more underestimation functions of f will be studied by applying the eigendecomposition approach to the D.C. forms provided by (i) and (ii) of Theorem 3. Specifically speaking, each of them can be rewritten as follows: T symmetric posi- tive semidefinite matrices.Hence, there exist two orthonormal matrices Ũ, Û ∈ ℝ n×n and two diagonal matrices D, D ∈ ℝ n×n such that Q 1 = Ũ D ŨT and Q 2 = Û D ÛT .The diagonal elements of D are the nonnegative eigenvalues λ1 , … , λn ∈ ℝ of Q 1 , while the orthonormal columns ũ1 , … , ũn ∈ ℝ n of Ũ are the corresponding eigenvectors of Q 1 ; in similar way, the diagonal elements of D are the nonnegative eigenvalues λ1 , … , λn ∈ ℝ of Q 2 , while the orthonormal columns û1 , … , ûn ∈ ℝ n of Û are the corresponding eigenvectors of Q 2 .As a consequence, since Q 1 and Q 2 have no negative eigenvalues, it results: 38 Page 18 of 32 for all i ∈ Θ + , and vi = √ λi ⋅ ûi for all ∀i ∈ Γ + .Hence, the following further D.C.
forms hold: In this light, the following branching variables are suggested by ( 4) and ( 5), respectively: so that For the sake of convenience, let Moreover, for all L ∶= ( L i ) i∈Θ + and U ∶= ( U i ) i∈Θ + such that L ≤ U and for all L ∶= ( L i ) i∈Γ + and U ∶= ( U i ) i∈Γ + such that L ≤ U , the following rectan- gles are introduced: Rectangles , and , result to be the smallest partitions Π 0 containing S with respect to branching variables and , respectively.
Theorem 6 Let f ∶ ℝ n → ℝ be expressed as in (i) and (ii) of Theorem 3.Then, the following convex quadratic functions are underestimation functions for f over S ∩ L , U and S ∩ L , U , respectively: and (5) 2 the maximum errors for Er 6 (x) and Er 7 (x) , respectively, obtained at ṽT i x = Proof The thesis follows in the same lines of Theorems 4 and 5. ◻

A particular case
In many applicative problems (see the forthcoming Sect.5), the linear multiplicative objective function f ∶ ℝ n → ℝ has the following particular structure: where z ∈ ℝ m , ∶= ( i ) i=1,…,p ∈ ℝ p , g ∶= (g i ) i=1,…,p ∈ ℝ p , q(z, , g) a linear or a convex quadratic term, and n ∶= m + 2p .In this light, it is worth studying the behavior of the underestimation functions proposed so far in the particular case of objective functions of type (6).

A computational experience
The solution method described and discussed in the previous sections has been implemented in a macOS 12.5.1 environment with an M1 Pro 10-core processor, MATLAB 2022a for coding and Gurobi 9.5.2 as solver for LP and QP problems.In this section, the results of some computational tests are presented, where the performances are compared with respect to the various underestimation functions previously studied.In this light, various instances have been randomly generated by using the "randi()" MATLAB function (integer numbers generated with uniform distribution).The average times spent to solve the instances (obtained with the "tic" and "toc" MATLAB commands), as well as the average number of iterations in procedure "Solve()" needed to solve them, are given as results of the computational tests.
The used tolerance parameters are relTol = 2 −35 and errTol = 2 −20 .For the sake of simplicity, in the instances generation we fixed the values a 0 = 0 and c 0i , d 0i = 0 for all i ∈ 1, … , p , we considered no equality constraints A eq x = b eq and a number of inequality constraints m = 2n .Moreover, the two following cases are taken into account in the instances generation:

Remark 6
The "nonnegative" case is aimed to study the behavior of the underestimation functions when both variables and branching variables are nonnegative, just like sometimes assumed in the literature (see, e.g., Zhou et al. 2015;Shen et al. 2020) and thus covering as a particular case the applicative problems described in Sect. 4.

A first comparison of all the underestimations
First of all, it is worth comparing all the introduced underestimation functions.
The standard value = 0.5 has been assumed for splitting all the underestimations.Starting from instances having n = 10 and p = 4 , the behaviors with p increased to p = 6 and with n increased to n = 20 are considered too.Moreover, both the "gen- eral" and the "nonnegative" cases are taken into account.The average times and average iterations are given as results of this first computational test and are summarized in the following tables.Six groups of instances (depending on n, p, and "general"/"nonnegative" cases) and 100 instances for each group are considered, with a grand total of 4200 problems solved.
The results provided in Tables 1 and 2 point out that: • In the "general" case the linear underestimation functions provide very bad tightness and hence very bad performance, while in the "nonnegative" case their performance has the same order of magnitude of the other underestimations; • Among the quadratic underestimation functions, u 4 is always much better than u 3 ; • Among the hybrid underestimation functions, u 7 is always much better than u 6 (this follows being u 6 derived from u 3 and being u 7 derived from u 4 ); • The most performing underestimation functions are always u 0 , u 4 , u 7 ; • In the "nonnegative" case, the performance of the underestimation functions are always better than the "general" case: this is due to the tightness of underestimation functions which results to be more effective in the "nonnegative" case than in the "general" one; • Among the linear and quadratic underestimation functions, increasing the value of p affects the performances much more than increasing the number of variables n; • Performances follow the number of branching variables of the various underestimation functions; in this light, recall that u 1 , u 2 and u 3 have 2p branching vari- ables, u 4 has p branching variables, while u 0 , u 6 and u 7 have |Λ − | , |Θ + | and |Γ + | branching variables, respectively.

A deep comparison of u
The previous subsection pointed out that the most performing underestimations are u 0 , u 4 and u 7 (and recall that, in the particular case of Sect. 4,these under- estimations coincide).The aim of this subsection is to focus on the behavior of these underestimations with respect to the parameter used to split the partitions.Assuming a number of variables n = 25 , instances for p = 4 , p = 7 and p = 10 are considered in both the "general" and the "nonnegative" cases.Values from 0 to 1 are tested.The average times and average iterations are given as results of this second computational test and are summarized in the following tables.Six groups of instances (depending on p, and "general"/"nonnegative" cases) and 100 instances for each group are considered, with a grand total of 19800 problems solved.Numbers in bold emphasize the best results (lower values) in terms of iterations or computational time.Notice that, at the best of our knowledge, no detailed studies have been published regarding to the impact of the splitting parameter in the behavior of the branch-and-bound method (for example, just = 0.25 is used in Gerard et al. (2017) and just = 0.8 is considered in Fampa et al. (2017)).Taking into account the results in Tables 3, 4, 5 and 6, it is worth noticing that: • underestimation functions u 4 and u 7 have similar performances when p < n ; in this light, notice that u 4 can be easily obtained while u 7 needs some eigenvec- tors to be computed; take into account also that in the case p > n underestima- tion u 7 has less branching variables than u 4 and hence is more performing; • u 0 has the best performances, but needs some eigenvectors to be computed; • the use of = 1.0 should be avoided since provides bad performances; • the greater is the value of p, the smaller is the value of providing the best performance; in this light, the parameters suggested in Fampa et al. (2017), Gerard et al. (2017) seem no useful; • performances in the "nonnegative" case are much better than the ones in the "general" case (underestimations' tightness results better in the "nonnegative" case than in the "general" one); • performances decrease exponentially with respect to the number of branching variables (and, hence, with respect to p); • as regards the particular class of problems described in Sect.4, in which u 0 , u 4 and u 7 coincide, the best choice is to use u 4 which needs no eigendecompositions.

A deep comparison of u
The aim of this subsection is to focus on the behavior of underestimations u 0 , u 4 and u 7 with respect to the number of variables n.Assuming a parameter p = 10 , instances for n = 5 , n = 10 and n = 20 are considered in both the "general" and the "nonnegative" cases.Values from 0 to 1 are tested.The average times and average iterations are given as results of this second computational test and are summarized in the following tables.Six groups of instances (depending on n, and "general"/"nonnegative" cases) and 100 instances for each group are considered, with a grand total of 19800 problems solved.Numbers in bold emphasize the best results (lower values) in terms of iterations or computational time.In other words, a detailed computational experience is provided to show the behavior of underestimations u 0 , u 4 and u 7 in the cases n < p , n = p and n > p .The results provided in Tables 7, 8, 9 and 10 point out that: • Underestimation u 0 is the most performing one (but needs eigenvalues and eigen- vectors to be computed); • Comparing u 4 and u 7 (recall that u 7 is derived from u 4 by means of an eigen- decomposition), performances are similar when n > p , u 7 is better than u 4 when n = p , while u 7 outperforms u 4 when n < p ; this behavior yields from the number of splitting variables which results to be smaller than or equal to min{n, p}; • Performances in the "nonnegative" case are much better than the ones in the "general" case; • The higher is the value of n the smaller the value of should be.

Overall comments
The main results of this computational experience are: • u 4 is the most performing underestimation function based on the structure of lin- ear multiplicative functions thus avoiding the numerical troubles of eigenvectors computing; • Linear underestimations, often used in the literature, are actually worse than quadratic underestimations; • The parameter has no value better than others, unlike the literature proposes; • In the particular applicative case described in Sect. 4 u 4 is the best choice and there is no need at all to use eigenvectors; • The "nonnegative" case results to have better performances than the "general" one, and this deserves to be deepened on in future researches.

Toward applications to bilevel problems
In Sect.3.4, a special linear multiplicative problem is considered in order to point out its behavior with respect to the underestimation functions previously introduced.At the same time, the study of this particular case is also motivated by the fact that the structure ( 6) is commonly present in several bilevel programming problems of leader-follower type (see, e.g., Dempe 2020 for a state of the art).Specifically speaking, leaderfollower type problems are hierarchical mathematical formulations with an upper and lower-level structure: the upper/leader-level is a suitable optimization problem partially constrained by a second (parametric) optimization problem as lower/follower-level.In many real-world applications formulated as bilevel programming problems, the upper/leader-level objective function is defined as in (6), that is, f (z, , g) = ∑ p i=1 i g i + q(z, , g) .By using the notation introduced in Sect. 3.4, if ∶= ( i ) i=1,…,p ∈ ℝ p represents the upper/leader variable (for instance, a vector of clearing prices), g ∶= (g i ) i=1,…,p ∈ ℝ p represents the lower/follower variable (for instance, a vector of certain quantities sold) and x ∶= (z, ) so that f (z, , g) ≡ f (x, g) , then the standard form of the resulting optmistic bilevel programming problem is the following: where and are the upper/leader and the lower/follower feasible sets, respectively.Under suitable assumptions (see, e.g., Dempe andZemkoho 2012, 2013) and under opportune constraints qualification conditions (see, e.g., Aussel and Svensson 2019; Dempe 2020), the optimistic bilevel programming problem (7) can be transformed into a (7) min x,g f (x, g) s.t.
x ∈ X g ∈ g ∈ K(x) ∶ h(x, g) = min g∈K(x) h(x, g) , single-level optimization problem1 , that is, a Mathematical Program with Equilibrium Constraints (MPEC) in which the lower/follower problem is reformulated by using opportune Karush-Kuhn-Tucker (KKT) optimality conditions.Notice that, in this case, the lower-level KKT conditions result to be: Then, the resulting KKT reformulation of the optimistic bilevel programming problem ( 7) is: Remark 7 The KKT reformulation allows to transform the optimistic bilevel programming problem ( 7) into a single-level program (9) at the cost of introducing the new variables in the formulation of the problem.In addition, the presence of complementarity constraints implies that the feasible region of the resulting problem is no more a polyhedron.
From a computational point of view, the linear multiplicative function f can be studied by using opportunely the results obtained in Sect.3.
Moreover, the complementarity constraints of the lower/follower problem in the KKT reformulation merged with the upper/leader-level, can be managed in various ways: • If the upper and lower level of (7) are linear multiplicative, then a classical procedure to solve bilevel programs is to consider a suitable penalization of T (x, g) that could be used to get a linear multiplicative single-level optimization problem (see Section 2.4.3 in Bard 1997, for example) at the cost of introducing the new penalty parameter in the formulation of the problem; • An outer approximation branch-and-bound method based on a feasible region relaxation where some of the products i l i l (x, g) = 0 are omitted and some of the variables i l or i l (x, g) are fixed to zero; • A quadratic indefinite penalty function M ∑ k l i l =1 i l i l (x, g) , with M >> 0 great enough, to be added to the objective function (again, a branch-and-bound may deserve); • Binary 0 − 1 variables and the so called big-M method (constraints i l i l (x, g) = 0 substituted with i l ≤ i l M and i l (x, g) ≤ (1 − i l )M , i l ∈ {0, 1} and M >> 0 (8) ∇ g h(x, g) + ∇ g (x, g) T = 0, ≥0, − (x, g)≥0, T (x, g) = 0. great enough), a branch-and-bound approach may be needed to manage the binary variables.
Remark 8 Complementarity constraints can be directly managed by means of a branch and bound approach where various subproblems are solved.Moreover, further approaches can be used in the case the available solvers are able to manage particular constraints.In this light, as i l ≥ 0 and − i l (x, g) ≥ 0 , then the following conditions (i)-(v) are equivalent: (i) i l i l (x, g) = 0; (ii) i l + i l (x, g) 2 = i l − i l (x, g) 2 ; (iii (iv) i l + i l (x, g) = max i l − i l (x, g);0 + max i l (x, g) − i l ;0 ; (v) i l , i l (x, g) is a "SOS1" set of variables (Special Ordered Set of type 1, set of variables where at most one variable may be nonzero, all others being at 0).
As a consequence, solvers able to manage quadratic indefinite constraints, absolute value constraints, "max" constraints or SOS constraints could be useful too.
In the context of this class of problems, although opportune numerical experiments could be of great interest to provide a performance comparison with works on the topic available in the literature (see Kleinert and Schmidt 2021, for example), it is beyond the scope of this work.Anyway, future investigations will be in this direction in relation to the study of suitable real-world problems.

Conclusions
In this paper, an extended computational experience regarding linear multiplicative problems is provided and a unifying framework to approach such problems with a branch-and-bound method is fully described.In this light, several underestimation functions are studied and various partitioning criteria are compared.In particular, it has been shown that the quadratic underestimation function u 4 has to be chosen to avoid eigenvectors-based approaches.As regards the splitting parameter , as higher the number of branching variables (or the number of variables) as smaller the value of should be.The proposed solution method results to be very efficient in the case of nonnegative variables and nonnegative branching variables.In this light, a particular case (useful in applications and in bilevel programming) has been theoretically studied in deep, pointing out also that u 4 is the underestimation to be used.The results obtained within the proposed unifying framework provide detailed comparisons and improvements with respect to the current literature.

Page 4 of 32 Definition 2
Let S ⊆ ℝ n be nonempty.Let f ∶ S → ℝ and Φ ∶ S → ℝ .Then, Φ is an underestimation function of f if:

•
A "general" case, where vectors a, , c i and d i have been randomly generated with components in the interval [−4, 4] , vectors l b and u b have been generated with components in the interval [−10, 10] , matrix A in has been generated with components in the interval [−10, 10] , b in has been generated in order to guarantee a feasible region different from the box [l b , u b ] and nonempty; • A "nonnegative" case, where vectors c i and d i have been randomly generated with components in the interval [0, 4], vector a has been generated with components in the interval [−4, 4] , vectors l b and u b have been generated with com- ponents in the interval [0, 15], matrix A in has been generated with components in the interval [−10, 10] , b in has been generated in order to guarantee a feasible region different from the box [l b , u b ] and nonempty.