Approximation of the Constant in a Markov-Type Inequality on a Simplex Using Meta-Heuristics

Abstract

Markov-type inequalities are often used in numerical solutions of differential equations, and their constants improve error bounds. In this paper, the upper approximation of the constant in a Markov-type inequality on a simplex is considered. To determine the constant, the minimal polynomial and pluripotential theories were employed. They include a complex equilibrium measure that solves the extreme problem by minimizing the energy integral. Consequently, examples of polynomials of the second degree are introduced. Then, a challenging bilevel optimization problem that uses the polynomials for the approximation was formulated. Finally, three popular meta-heuristics were applied to the problem, and their results were investigated.

1. Introduction

Markov-type polynomial inequalities (or inverse inequalities) [1,2,3,4,5,6], similarly to Bernstein-type inequalities [1,5,6,7,8,9], are often found in many areas of applied mathematics, including popular numerical solutions of differential equations. Proper estimates of optimal constants in both types of inequalities can help to improve the bounds of numerical errors. They are often used in the error analysis of variational techniques, including the finite element method or the discontinuous Galerkin method used for solving partial differential equations (PDEs) [10]. Moreover, Markov’s inequality can be found in constructions of polynomial meshes, as it is a representative sampling method, and is used as a tool for studying the uniform convergence of discrete least-squares polynomial approximations or for spectral methods for the solutions of PDEs [11,12,13]. Generally, finding exact values of optimal constants in a given compact set E in ${\mathbb{R}}^{\mathbb{N}}$ is considered particularly challenging. In this paper, an approach to determining the upper estimation of the constant in a Markov-type inequality on a simplex is proposed. The method aims to provide a greater value than that estimated by L. Bialas-Ciez and P. Goetgheluck [4].

For an arbitrary algebraic polynomial P of N variables and degree not greater than n, Markov’s inequality is expressed as

$${\parallel \left|gradP\right|\parallel }_E\le C{n}^2{\parallel P\parallel }_E,$$

(1)

where $gradP:=\left(\frac{\partial P}{\partial x_1},\dots ,\frac{\partial P}{\partial x_N}\right),\left|gradP\right|:={\left({\sum }_{j=1}^N\left(\frac{\partial P}{\partial x_j}\right)\right)}^{1/2},$ ${\parallel P\parallel }_E=sup\left\{\right|P\left(x\right)|:x\in E\},$ $C>0$, E is a convex compact subset of ${\mathbb{K}}^N$ (with non-void interior), and $\mathbb{K}\in \mathbb{R}.$ The exponent 2 in (1) is optimal for a convex compact set E (with non-void interior) [2,4].

In 1974, D. Wilhelmsen [14] got an upper estimate for the constant $C<\frac{4}{{\omega }_E}$ in Markov’s inequality (1), where ${\omega }_E$ is the width of E. The width of $E\in {\mathbb{R}}^{\mathbb{N}}$ is defined as ${\omega }_E:=inf\{dist(P{l}_1,P{l}_2)\mid (P{l}_1,P{l}_2)\in G\}$ [4], where G is the set of all possible couples of parallel planes $P{l}_1$ and $P{l}_2$ such that E lies between $P{l}_1$ and $P{l}_2.$ Furthermore, for any convex compact set $E\subset {\mathbb{R}}^N$, the author proposed the hypothesis that $C=\frac{2}{{\omega }_E}.$ M. Baran [3,7] and Y. Sarantopoulos [15] independently proved that Wilhelmsen’s conjecture is true for centrally symmetric convex compact sets. In turn, L. Bialas-Ciez and P. Goetgheluck [4] showed that the hypothesis is false in the simplex case and proved that the constant $C>2.842534.$

The best constant in Markov’s inequality (1) for a convex set was also studied by D. Nadzhmiddinow and Y. N. Subbotin [6] and by A. Kroó and S. Révész [16]. It is worth mentioning that for any triangles in ${\mathbb{R}}^2$, A. Kroó and S. Révész [16] generalized the result for a simplex given in [4], considering an admissible upper bound for $C\le \frac{\sqrt{10}}{{\omega }_E}$, and received an upper estimate value of constant C in (1): $C\le \frac{s\left(\mathsf{\Delta }\right)}{w\left(E\right)},$ where $s\left(\mathsf{\Delta }\right):=\frac{2}{a}\sqrt{a^2+b^2+2abcos{\gamma }^{\prime }},$ with $\mathsf{\Delta }$ being a triangle with the sides $a,b,c$ and the angles ${\alpha }^{\prime },{\beta }^{\prime },{\gamma }^{\prime }$, respectively, such that $c\le b\le a$ and ${\gamma }^{\prime }\le {\beta }^{\prime }\le {\alpha }^{\prime }.$

In the recent literature, many works have considered Markov inequalities. For example, S. Ozisik et al. [10] determined the constants in multivariate Markov inequalities on an interval, a triangle, and a tetrahedron with the $L_2$-norm. They derived explicit expressions for the constants on all above-mentioned simplexes using orthonormal polynomials. It is worth noticing that exact values of the constants are the key for the correct derivation of both a priori and a posteriori error estimates in adaptive computations. F. Piazzon and M. Vianello [11] used the approximation theory notions of a polynomial mesh and the Dubiner distance in a compact set to determine error estimates for the total degree of polynomial optimization on Chebyshev grids of the hypercube. Additionally, F. Piazzon and M. Vianello [12] constructed norming meshes for polynomial optimization by using a classical Markov inequality on the general convex body in ${\mathbb{R}}^{\mathbb{N}}$. A. Sommariva and M. Vianello [13] used discrete trigonometric norming inequalities on subintervals of the period. They constructed norming meshes with optimal cardinality growth for algebraic polynomials on the sections of a sphere, a ball, and a torus. Moreover, A. Kroó [17] studied the so-called optimal polynomial meshes and O. Davydov [18] presented error bounds for the approximation in the space of multivariate piecewise polynomials admitting stable local bases. In particular, the bounds apply to all spaces of smooth finite elements, and they are used in the finite element method for general fully nonlinear elliptic differential equations of the second order. L. Angermann and Ch. Heneke [19] proved error estimates for the Lagrange interpolation and for discrete $L_2$-projection, which were optimal on the elements and almost optimal on the element edges. In another study, M. Bun and J. Thaler [5] reproved several Markov-type inequalities from approximation theory by constructing explicit dual solutions to natural linear programs. Overall, these inequalities are a basis for proofs of many approximated solutions that can be found in the literature.

In this paper, an approach to determining an approximation of the constant in Markov’s inequality using introduced minimal polynomials is proposed. Then, the polynomials are used in the formulated bilevel optimization problem. Since many optimization techniques could be employed to find the constant, three popular meta-heuristics are applied and their results are investigated.

The rest of this paper is structured as follows. In Section 2, minimal polynomials are introduced. Then, in Section 3, their application in the approximation of the constant in a Markov-type inequality is considered. In Section 4, in turn, the approximation is defined as an optimization problem and the problem is solved using three meta-heuristics. Finally, in Section 5, the study is summarized and directions for future work are indicated.

2. Minimal Polynomials

The monic polynomial is defined as: $P\left(x\right):=x^{\alpha }+{\sum }_{\beta \prec \alpha }{a}_{\beta }{x}^{\beta },$ for each $x\in {\mathbb{C}}^{\mathbb{N}}$, and with $a_{\beta }\in \mathbb{C}$ [20,21,22] (the coefficient of its leading monomial is equal 1), where $\alpha =({\alpha }_1,\dots ,{\alpha }_N)\in {\mathbb{N}}^N,$ $\left|\alpha \right|:={\alpha }_1+\cdots +{\alpha }_N$ is its length and $x^{\alpha }:=x_1^{{\alpha }_1}{x}_2^{{\alpha }_2}\cdots x_N^{{\alpha }_N};$≺ indicates the graded lexicographic order, defined as: If $\left|\alpha \right|<|{\alpha }^{\prime }|$, then $\alpha \prec {\alpha }^{\prime }$ $(x^{\alpha }\prec x^{{\alpha }^{\prime }})$, while if $\left|\alpha \right| = |{\alpha }^{\prime }|$ when $\alpha \prec {\alpha }^{\prime }(x^{\alpha }\prec x^{{\alpha }^{\prime }})$, the first (starting from the left) non-zero coefficient of ${\alpha }^{\prime }-\alpha$ is positive [21,22]. Polynomial $P\left(x\right)$ of the smallest deviation from zero on E with respect to the supremum norm is called the minimal polynomial (Chebyshev polynomial) for E [22,23,24,25].

Let E be a compact set in a complex plane that contains at least $\left|\alpha \right|+1$ points. Based on the classical theorem of Tonelli [26], among the monic polynomials of degree $\left|\alpha \right|$, there exists only one polynomial $T_{\left|\alpha \right|}$ that minimizes the supremum norm on E. The $T_{\left|\alpha \right|}$ is called the Chebyshev polynomial of order $\left|\alpha \right|.$ It is defined as $\parallel T_{\left|\alpha \right|}{\parallel }_E=inf\left\{\parallel x^{\left|\alpha \right|}+{\sum }_{i=0}^{\left|\alpha \right|-1}{a}_i{x}^i{\parallel }_E\right\}$, where the infimum is taken with respect to $(a_0,a_1,\dots ,a_{\left|\alpha \right|-1})\in {\mathbb{C}}^{\left|\alpha \right|}.$ In other words, $T_{\left|\alpha \right|}-x^{\left|\alpha \right|}$ is the polynomial of best approximation to $x^{\left|\alpha \right|}$ on $E.$ This polynomial is called the minimal (monic) polynomial of degree $\left|\alpha \right|.$

The complex equilibrium measure for a standard simplex in ${\mathbb{R}}^N$ is calculated for $E=\{x\in {\mathbb{R}}^N:x_1,\dots ,x_N\ge 0,x_1+\cdots +x_N\le 1\}=S_N$, which is a simplex in ${\mathbb{R}}^N$ [8]. Then, the plurisubharmonic extremal function on E, $u_E\left(z\right)=logh\left(\right|z_1|+\cdots +|z_N|+|z_1+\cdots +z_N-1\left|\right),$ is used to calculate the (complex) equilibrium measure ${\lambda }_{S_N}=N!vol\left(B_N\right){(x_1·\cdots ·x_N)}^{-1/2}{(1-x_1-\cdots -x_N)}^{-1/2}dx$, where $vol\left(B_N\right)$ is the usual volume of the closed-unit Euclidean ball in ${\mathbb{R}}^{\mathbb{N}}$.

3. Estimation of the Constant $\mathit{C}$ in a Markov-Type Inequality

Let $S_2$ be a simplex in ${\mathbb{R}}^2.$ Consequently,

$$S_2=\{(x_1,x_2)\in {\mathbb{R}}^2:x_1,x_2\ge 0,x_1+x_2\le 1\}.$$

The normalized equilibrium measure $S_2$, as a subset of ${\mathbb{C}}^2$, is a two-dimensional Lebesgue measure with the weight $A(x_1,x_2):=\frac{1}{\sqrt{x_1}\sqrt{x_2}\sqrt{1-x_1-x_2}}.$ For continuous functions $f,g$ on $S_2$, the scalar product is defined as

$$\left(f,g\right)=\underset{S_2}{\int \int }f(x_1,x_2)g(x_1,x_2)A(x_1,x_2)d{x}_{1}d{x}_2.$$

Then, a polynomial $Q(x_1,x_2)=x_1+x_2-1/2$ is defined. It is the minimal polynomial for the simplex, and $P(x_1,x_2)=a_0{x}_1^2+a_1{x}_2^2+a_2{x}_1{x}_2+a_3{x}_1+a_4{x}_2+a_5,$ for some real numbers $a_0$, $a_1$, $a_2$, $a_3$, $a_4$, and $a_5.$ Thus,

$$\begin{array}{c}(P,Q)=\underset{S_2}{{\displaystyle \int \int }}P(x_1,x_2)Q(x_1,x_2)A(x_1,x_2)d{x}_{1}d{x}_2=a_0{x}_{3,0}+a_1{x}_{1,2}+a_2{x}_{2,1}+a_3{x}_{2,0}+a_4{x}_{1,1}\\ +a_5{x}_{1,0}+a_0{x}_{2,1}+a_1{x}_{0,3}+a_2{x}_{1,2}+a_3{x}_{1,1}+a_4{x}_{0,2}+a_5{x}_{0,1}-a_0{x}_{2,0}/2+a_1{x}_{0,2}/2\\ -a_2{x}_{1,1}/2-a_3{x}_{1,0}/2-a_4{x}_{0,1}/2-a_5{x}_{0,0}/2,\end{array}$$

where $x_{i,j}$ for $i,j\in \{0,1,2,3\}$ are symbols given in

Table 1. Symbols and their values.

SYMBOLS	${\mathit{x}}_{0,0}$	${\mathit{x}}_{1,0}$	${\mathit{x}}_{0,1}$	${\mathit{x}}_{1,1}$	${\mathit{x}}_{1,2}$	${\mathit{x}}_{2,1}$	${\mathit{x}}_{2,0}$	${\mathit{x}}_{0,2}$	${\mathit{x}}_{3,0}$	${\mathit{x}}_{0,3}$
VALUES	$2\pi$	$2\pi /3$	$2\pi /3$	$2\pi /15$	$2\pi /35$	$2\pi /35$	$2\pi /5$	$2\pi /5$	$2\pi /7$	$2\pi /7$

.

We want the polynomial P to be orthogonal to Q so that $(P,Q)=0.$ Substituting the values given in Table 1 into the expression above, $a_0/7+a_1/7+a_2/21+a_3/5+a_4/5+a_5/3=0$ is obtained. Finally, $a_5=-3a_0/7-3a_1/7-a_2/7-3a_3/5-3a_4/5.$

Here, a counterexample to the Wilhelmsen hypothesis can be formulated, which is different from that discussed by L. Bialas-Ciez and P. Goetgheluck [4], to support the claim that the approach leads to a better approximation. Counterexample: Consider the following polynomial $P(x_1,x_2)$ $=a_0{x}_1^2+a_1{x}_2^2+a_2{x}_1{x}_2+a_3{x}_1+a_4{x}_2+a_5$, where $a_5=-3a_0/7-3a_1/7-a_2/7-3a_3/5-3a_5/5.$ Putting $a_0=6,a_1=3/4,a_2=33/4,a_3=-6,a_4=-57/28,a_5=3/4$, $P(x_1,x_2)=6x_1^2+0.75x_2^2+8.25x_1{x}^2-6x_1-2.04x_2+0.75$ is obtained. It is easy to calculate that ${∥P∥}_{S_2}$ is reached for point $\left(13/14,1/14\right)$, and its value is $297/392.$ Similarly, ${∥gradP(0,1)∥}_{S_2}=\sqrt{58500}/28.$ Hence, this simple and far from optimal approximation value is greater than those of L. Bialas-Ciez and P. Goetgheluck [4]: ${\parallel gradP\parallel }_{S_2}/{(4\parallel P\parallel }_{S_2})=2.850293143.$

4. Optimization Problem in a Markov-Type Inequality

The estimation of the constant C in Markov’s inequality can be improved by considering a minimal polynomial and an equilibrium measure. The exemplary calculations of minimal polynomials and their corresponding Markov inequalities on the simplex are shown in

Table 2. Constants C in Markov inequalities on the simplex for exemplary polynomials.

Example	Minimal Polynomial $\mathit{P}(\mathit{x},\mathit{y})$	${\parallel \mathit{P}\parallel }_{{\mathit{S}}_2}$	${\parallel \mathit{gradP}\parallel }_{{\mathit{S}}_2}$	${\parallel \mathit{gradP}\parallel }_{{\mathit{S}}_2}/4{\parallel \mathit{P}\parallel }_{{\mathit{S}}_2}$
1	$y^2-y+\frac{1}{9}+xy-\frac{4}{9}{x}^2+\frac{2}{9}x$	5/36	$\sqrt{202}$/9	C = 2.842534081
2	$y^2-y+\frac{1}{8}+\frac{1}{8}{x}^2+\frac{33}{24}xy-\frac{19}{56}$	99/784	10$\sqrt{585}$/168	C = 2.850293143
3	$y^2-y+\frac{1}{8}+\frac{61}{40}xy+\frac{1}{4}{x}^2-\frac{1}{2}x$	221/1760	$\sqrt{3281}$/40	C = 2.851041419
4	$y^2-y+\frac{1}{8}+\frac{3}{2}xy+\frac{1}{4}{x}^2-\frac{1}{2}x$	1/8	$\sqrt{2}$	C = 2.828427125

. In the first example of the polynomial, the obtained constant is equal to the value reported by L. Bialas-Ciez and P. Goetgheluck [4]. In the remaining two polynomials, the constants are improved. However, a slight change in the coefficients decreases the constant, as is shown in the fourth example. Therefore, in this paper, the polynomial is proposed as a result of a nontrivial optimization problem in which an improvement of the estimate of the constant C is considered.

In this paper, the constant in Markov’s inequality on the simplex for the proposed minimal polynomials is estimated. Assuming that

$${\parallel gradP(x,u)\parallel }_{S_2}/{(4\parallel P(x,u)\parallel }_{S_2})=C(x,u),$$

(2)

where P is one of the following minimal polynomials:

$$P_1=x_2^2-a_0{x}_2+a_1+a_2{x}_1{x}_2+a_3{x}_1^2+a_4{x}_1,$$

(3)

$$P_2=x_1^2-a_0{x}_1+a_1+a_2{x}_1{x}_2+a_3{x}_2^2+a_4{x}_2,$$

(4)

$$P_3=x_1^2-x_1+1/8+a_0{x}_1{x}_2+a_1{x}_2^2+a_2{x}_2,$$

(5)

$$P_4=x_2^2-x_2+1/8+a_0{x}_1{x}_2+a_1{x}_1^2+a_2{x}_1,$$

(6)

and $u=[a_0,a_1,a_2,a_3,a_4]$, in order to determine the greatest value of the function $C(x,u),$ an optimization problem is formulated.

R. M. Aron and M. Klimek [27] gave explicit formulas for supremum norms: ${\parallel (a,b,c)\left|\right|}_{\mathbb{R}}$ $=sup\left\{\right|a{x}^2+bx+c|:x\in [-1,1]\}$ and ${\parallel (a,b,c)\left|\right|}_{\mathbb{C}}=sup\left\{\right|a{z}^2+bz+c|:z\in \mathbb{C},|z|\le 1\},$ where $(a,b,c)\in {\mathbb{R}}^3.$

The vectors $u=[a_0,a_1,a_2,a_3,a_4]$ and $x=[x_1,x_2]$ are regarded as decision variables in the optimization problem. Consequently, $\widehat{u}$ and $\widehat{x}$ denote the vectors that maximize the $C(x,u)$. Considering the relationship between $P_1$ and $P_2,$ or $P_3$ and $P_4,$ only the coefficients that maximize their corresponding $C_1(x,u)$ and $C_3(x,u)$ need to be found, i.e., they are then applied to the remaining polynomials. For polynomials $P_3$ and $P_4,$ $a_3=a_4=0.$ The calculation of Equation (2) requires finding the extrema of the polynomial P inside the simplex limited by values $(x_1,x_2)\in \{(0,0),(0,1),(1,0)\}$. Therefore, the following bilevel optimization problem is solved. Here, one optimization task is nested within the other, i.e., the decision variables of the upper problem are the parameters in the lower one, and the lower problem is a constraint to the upper problem, limiting the feasible u values to those that are lower-level optimal [28].

The problem is given by

$$\begin{array}{ccc}\hfill & \underset{u\in R^5,x\in [0,1]}{\max }\hfill & \hfill C(x,u)\\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill \widehat{x}=\underset{x\in [0,1]}{\mathrm{argmax}}\{C(x,u):x_1+x_2-1\le 0\},\end{array}$$

(7)

where $argmax$ denotes the arguments of the function that maximize its value.

Finally, the values of $\widehat{x}$ and $\widehat{u}$ are found as the solutions to the optimization problem. In order to determine the greatest value of C, the nonlinearity of the considered objective function, as well as the existence of multiple local optima, must be taken into account. Here, different values of u can lead to the same value of C, or close values of u can result in a large difference between the corresponding values of $C.$ This can be seen in Figure 1, in which only a small part of the search space is visualized. For the visualization, C is calculated, keeping one coefficient constant and selecting the remaining one from the range $[-1,1]$ with the step $0.1$. The coefficients are selected from the set $\{-0.1,0,0.1\}.$ The difference between the values of the constant is purposely small to show the changes of values of $C.$ Each subplot in Figure 1 contains a marker to additionally highlight such changes.

To solve the defined optimization problem, a suitable method should be determined. Therefore, in this paper, a hybrid approach is proposed. It is justified by the need to find the extrema of the polynomial P inside the simplex for the computation of C. Such extrema are found using the Interior Point Algorithm (IPT) [29,30], i.e., $\widehat{x}$ is found for an examined u from the upper level. The upper-level optimization is more difficult. Since the considered problem has not been formulated and solved in the literature, it cannot be determined in advance which method should be used to determine the solution. Therefore, taking into account possible time-consuming computations, three popular meta-heuristics are used and compared: the Particle Swarm Optimization method (PSO) [31,32,33], Genetic Algorithm (GA) [34,35], and Simulated Annealing technique (SA) [36,37]. The application of population-based algorithms for solving different bilevel problems is popular (e.g., [28]).

The PSO is a global evolutionary method based on swarm intelligence. It finds the optimal solution by exploring the search space using a set of particles. Each particle represents a candidate solution to the problem. The PSO mimics the movement of a bird flock or a fish school; hence, the dynamics of the set of particles in the solution space are iteratively governed by the following equations [31]:

$$\begin{array}{c}\hfill p_i\left(n\right)=p_i(n-1)+v_i\left(n\right)\\ \hfill v_i\left(n\right)=w{v}_i(n-1)+c_{1}rn{d}_1(p{b}_i-p_i(n-1))+c_{2}rn{d}_2(pB-p_i(n-1)),\end{array}$$

where $p_i\left(n\right)$ and $v_i\left(n\right)$ are the position and the velocity of the i-th particle at time $n,n=1,2,\dots ,N$; $p{b}_i$ is the best position of the i-th particle, while $pB$ is the best position of the swarm at n. w is the inertia weight, and $c_1$ and $c_2$ are coefficients constraining the parts of the equation that simulate cognitive and social factors of the swarm [31,32,33]); $rn{d}_1$ and $rn{d}_2$ are random numbers drawn from the range $[0,1]$. The termination criterion of the algorithm is the number of iterations N or the violation of the minimal acceptable difference in the objective function C between consecutive iterations. Both criteria are commonly used in most heuristic optimization algorithms. Finally, the PSO returns $\widehat{u}=pB\left(N\right)$.

The GA uses a population of individuals (called chromosomes), representing candidate solutions [34,35]. In the GA, at each generation (iteration), a population of individuals is assessed to select the best solutions, which are then used to create the population in the next generation. The assessment is based on the values of the objective function, while new individuals are obtained using crossover and mutation operators. The crossover combines the information carried by two selected parents in one solution by mixing values of their vector representations. The mutation, in turn, introduces random changes to the individual. After the predefined number of generations, the algorithm returns the best individual $\widehat{u}$.

The simulated annealing, unlike the population-based PSO and GA, uses much fewer resources and requires fewer objective function computations, since it processes only one solution at each iteration [36,37]. The SA mimics the process of heating a metal and slowly lowering its temperature. The algorithm starts from a random solution, which is iteratively modified. At each iteration, reflecting a temperature drop, the algorithm modifies the current solution and compares the values of the objective function before and after the modification. If a new solution is better or the difference between their values is within a predefined range, the modification is accepted. The acceptance of a slightly worse solution helps the method to escape local minima [36,37].

All of these methods are commonly used for solving nonlinear optimization problems with local optima. However, their suitability for a given problem is mostly determined experimentally. In this paper, they are used to determine $\widehat{u}$. However, in order to calculate C for any u, the lower-level optimum $\widehat{x}$ must be solved each time. In this paper, it is obtained using the IPT algorithm [29,30]. Finally, the considered problem is solved using the hybrid approach [28]. In the IPT algorithm, an optimization problem is defined as

$$\begin{array}{ccc}\hfill & \underset{\mathit{x}}{\min }f\left(\mathit{x}\right)\hfill & \hfill \\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill h\left(\mathit{x}\right)=0andg\left(\mathit{x}\right)\le 0,\end{array}$$

(8)

where $f\left(\mathit{x}\right)$ is objective function to be minimized and h and g are equality and inequality constraints, respectively. The solved approximate problem is defined as

$$\begin{array}{ccc}\hfill & \underset{\mathit{x},\mathit{s}}{\min }{f}_{\mu }(\mathit{x},\mathit{s})\hfill & \hfill =\underset{\mathit{x},\mathit{s}}{\min }f\left(\mathit{x}\right)-\mu \sum _{i}ln\left(s_i\right)\\ \hfill & \mathrm{subject}\mathrm{to}\hfill & \hfill h\left(\mathit{x}\right)=0andg\left(\mathit{x}\right)+s=0,\end{array}$$

(9)

where $s_i$ are slack variables and $\mu >0$ is a barrier parameter. As $\mu$ converges to zero, the sequence of approximate solutions converges to the minimum of f [29,30].

In the numerical experiments, the solution processed by heuristic algorithms was represented by a five-dimensional vector of floating-point values, since $u=[a_0,a_1,a_2,a_3,a_4]$, and the IPT used a two-dimensional vector $x=[x_1,x_2]$. To ensure a fair comparison of meta-heuristics, each technique was run 50 times, and a maximum of 25,000 objective function calls were performed in each run. The number of function calls was achieved by selecting a large-enough population size, aiming at the best performance of the PSO and GA. The implementations of the algorithms that can be found in the Matlab Optimization Toolbox are used. Algorithms were run with their default parameters: the PSO (w = [01, 1.1], $c_1=1.4$, $c_2=1.4$), the GA (EliteCount = 0.05 of the population size, CrossoverFraction = 0.8, the scattered crossover, Gaussian mutation, and stochastic uniform selection), and the SA (InitialTemperature = 100, temperature is equal to InitialTemperature * 0.95${}^{iteration}$). Note that the values in u were not bounded (i.e., $u\in {\mathbb{R}}^5$), and the GA and PSO, unlike the SA, did not require the specification of the starting point of the optimization, as the initial populations of the solutions were selected randomly. For this reason, the SA started from $u=[0,0,0,0,0]$. As demonstrated in

Table 3. Optimization results obtained by the compared algorithms. The best values of C obtained in the first experiment are written in bold.

Algorithm	${\mathit{C}}_1,{\mathit{C}}_2$	${\mathit{C}}_3,{\mathit{C}}_4$
GA	2.861879141	2.852924202
SA	2.917638439	2.854084285
PSO	2.925401560	2.854206134

, the PSO was able to determine the best C for each polynomial. The algorithms explored the solution space in different ways, obtaining different results.

To further show the difficulty of the optimization task for the used algorithms, Figure 2 presents the convergence of the algorithms towards the optimum. Here, the progress of the set of solutions at each calculation of the objective function is reported. As presented, the approaches with the PSO and GA quickly converge, and most of the time is spent on the exploration of the solution space near the optimum. The SA, unlike the population-based methods, processes one solution at each iteration, which influences its convergence rate. The results reveal that the PSO is better suited to the problem than the GA, since it can determine a better C (see Table 3), and it finds it with fewer function calls. Interestingly, the set of solutions in the first function call is more diverse for the PSO, which may contribute to its superior performance, as in the remaining computation time, it is more focused on escaping the local optima. Similarly, it can be seen that the solutions emerging in the run of the SA exhibit large variability. Note that the number of runs (50) and the maximal number of function calls in these experiments (25,000) are purposely greater than could be concluded from the convergence plots to give the methods enough resources to find better solutions. However, the methods often terminate the computations if the change in the best function value is too small.

This experiment indicates that the proposed approach with the PSO can provide the best C. Therefore, in the second experiment, the parameters of this algorithm were set to conduct an even better search of the solution space. By increasing the number of particles and generations, the maximal number of function calls was set to 25 million, and the method was run 100 times. Note that the estimated time to finish one run of the PSO with 25 million function calls in Matlab on the i7-6700 K 4 GHz with 64 GB RAM CPU is 335 h. The calculations were run in a parallel manner, and the methods were allowed to terminate the optimization if a function value change was small. They are presented in

Table 4. Values of C obtained in the second experiment and their corresponding coefficients.

$C_1,C_2$	2.925415050
$a_0$	−2765.096993
$a_1$	−1383.048545
$a_2$	−15148.04101
$a_3$	−10633.59420
$a_4$	10846.85282
$C_3,C_4$	2.854206170
$a_0$	1.169292858
$a_1$	−0.121651344
$a_2$	−0.129695921

. The second experiment improved the already obtained values for $C_1$ (or $C_2$) and $C_3$ (or $C_4$) at the fifth and eighth decimal places, respectively (see Table 3). The run in which the best $C_1$ was obtained required 720 thousand function calls, and $C_3$ was found in a run after the objective function was calculated 1570 thousand times.

5. Conclusions

In this paper, it was demonstrated that the best constant C in Markov’s inequality for the minimal polynomials on the simplex can be determined as a solution to a nested optimization problem. First, the examples of the polynomials that allow estimation of the constant were proposed. Then, the optimization problem was defined and its solution was investigated. Since it could not be determined in advance which heuristic optimization algorithm should be used, three popular meta-heuristics were adopted and compared. Finally, the optimal value of the constant was obtained with a hybrid of the ITP and PSO.

Future work will focus on the application of the polynomials of the higher degrees to the considered problem. In addition, the usage of the introduced bilevel optimization problem as an evaluation tool for the comparison of the performance of meta-heuristics seems to be a promising future direction.