ITSO: A novel Inverse Transform Sampling-based Optimization algorithm for stochastic search

Despite the numerous research works and industrial applications of Artificial Intelligence algorithms, they have been criticized about lacking a solid theoretical background [1]. The empirical results demonstrate impressive performance, however, their theoretical foundation and analysis are often vague [2]. Machine Learning (ML) models, frequently aim at identifying an optimal solution [3, 4], which is computationally hard and attempts to explain the procedure often based on the evaluation of the function’s Gradients [5, 6]. Accordingly, the identification of whether a stochastic algorithm will work or not remains an open question for many research and real-world applications [7, 8, 9, 10]. A vast number of optimization algorithms have been developed and applied for the solution of the corresponding problems [11, 12, 13, 14], as well as mathematical proofs regarding their algorithmic convergence [15, 16, 17], however, their basic formulation often stems from nature-inspired procedures [18, 19] and not solid mathematical frameworks. Accordingly, a vast number of research works have been published in order to investigate the performance of black-box algorithms [20, 21, 22].

In its elementary form, the purpose of efficient optimization algorithms is to find the argument yielding the minimum value of a black-box function $f(x)$, defined on a set $A$, $f:A\to\mathbb{R}^{n}$. Accordingly, the inverse problem of maximization is the minimization of the negated function $-f(x)$, while problems with multiple objective functions often utilize single function optimization algorithms to attain the best possible solution. $A$ is assumed a compact subset of the Euclidean space ${{\mathbb{R}}^{n}}$, where $n$ is the number of dimensions of the set $A$, however, the proposed method applies similarly to discrete and continuous topological spaces. The unknown, black-box function $f$, returns values for the given input $x_{ij}={{\mathbf{x}}_{i}}=({{x}_{i1}},{{x}_{i2}},\cdots,{{x}_{in}})$ at each computational discrete time step $i=\left\{1,2,\ldots,{{f}_{e}}\right\}$, where ${{f}_{e}}$ is the number of maximum function evaluations. The sought solution is a vector ${{\mathbf{x}}_{min}}\in A$, such that $f({{\mathbf{x}}_{min}})\leq f(\mathbf{x}),\forall\mathbf{x}\in A$, which may be written by

The purpose of this work is to provide a rigor context for the optimization problem of a black-box function, by adhering to the Probability Theory, aiming at identifying the best possible solution $\mathbf{x}_{min}$, within the given iterations $f_{e}$, during the execution of the algorithm.

The rest of the paper is organized as follows: In Section 2, the proposed Inverse Transform Sampling Optimizer (ITSO) is presented. Additionally, the same Section provides in details the theoretical formulation of the algorithm as well as some programming and implementation techniques. Illustrative examples of the optimization history are also comprised. In Section 3 the theoretical proof of convergence is provided, as well as Lemma 1, deriving that the suggested optimization framework is the fastest possible. The numerical experiments are divided into three groups. Subsection 5 is about the comparison with 13 nonlinear loss functions, 17 optimization methods, for 10 and 20 dimensions of search space, and 5000 and 10000 iterations per dimension. Section 4, briefly presents the programming techniques that were investigated, in order to implement the proposed method into a computer code. Finally, the conclusions are drawn in Section 6.

Let the probability distribution of the optimal values of $f$, be considered as the product of some monotonically decreasing kernel $k$ over $f$ at some time-step $i$ and hence given input ${{\mathbf{x}}_{i}}$. This assumption stands in the foundation of the method and can be considered as rational, instead of an arbitrary selection strategy, inspired by natural or other phenomena. It is a straightforward application of the Probability theory to the problem. The kernel may be considered as parametric, concerning time $i$. The selection of the kernel $k$ should satisfy the condition that its limit is the Dirac delta function centered at ${\operatorname{arg\,min}}\,f$,

where $\delta_{m}$ denotes the Dirac function centered at ${\operatorname{arg\,min}}\,f$. Although the selection of the kernel $k$ is ambiguous, it can be chosen among a variety of functions satisfying Equation 2, such as the Gaussian:

where $g(i)$ is a time increasing pattern. The function $g(i)$ controls the shape of the kernel $k$, approximating numerically Equation 2. Additionally, we may use:

or any other non-negative Lebesgue-integrable function, which reverses the order of the given set of all ${{f}_{{\hat{i}}}}$, where $\hat{i}$ is the permutation of the indices $1,2,...,i$, such that the sequence of ${{f}_{{\hat{i}}}}$ being monotonically strictly decreasing. The duplicate values of ${{f}_{{\hat{i}}}}$ should be extracted to avoid numerical instabilities. These duplicates often appeared in the empirical calculations, especially when the algorithm was close to a local or global stationary point. A variety of kernel functions $k$ were investigated, and the results weren’t affected significantly, even when distorting $k(f)$ with some random noise ${X}^{\prime}\sim\mathcal{U}(a,b)\forall k({{f}_{{\hat{i}}}})\in(a,b)$.

Accordingly, for each dimension $j$ of the vector space $A$, the marginal probability density function ${P}_{{{X}_{j}}}$ is obtained numerically from the kernel function:

${{P}_{{{X}_{j}}}}({{x}_{ij}})dx$ is the probability that ${\operatorname{arg\,min}}\,f$ falls within the infinitesimal interval $[{{x}_{ij}}-dx/2,{{x}_{ij}}+dx/2]$. The corresponding cumulative distribution function (CDF) ${{F}_{{{X}_{j}}}}$, can be calculated by:

where $l{{b}_{j}}$ is the lower bound of the $j^{th}$ dimension of the vector space $A$ and can be numerically evaluated by some numerical integration rule, such as the Riemann sum $S_{j}=\sum\limits_{i=1}^{n}{{{P}_{{{X}_{j}}}}}(x_{ij}^{*})\Delta{{x}_{ij}}$, with ${{{P}_{{{X}_{j}}}}}(x_{ij}^{*})$ computed by the application of the kernel $k$ on some ${x}_{ij}$, such that $f(x_{ij}^{*})=(f({{x}_{ij}})+f({{x}_{i-1,j}}))/2$, or another approximation scheme. In the following pseudo-code 1, the algorithmic implementation of the Inverse Transform Sampling method is demonstrated. The symbols are also noted in the Nomenclature section.

A graphical demonstration of the evolution of the Probability Density, as well as the corresponding Cumulative Distribution Functions, is presented in Figure 1, for $f(x)=(x-5)^{2}$, which is function with one extreme value and in Figure 2, for $f(x)=sin(x+0.7)+0.01*(x+0.7)^{2}$, which has many extrema. Interestingly, as the function evaluations increase, the CDF, is characterised by sharp slopes, which are positioned in regions where the PDF exhibits high values, and hence function $f(x)$ attends its lows. Figures 1 and 2, offer an intuitive representation of the procedure for finding the minimum of the function, within the suggested framework.

During the optimization process, the algorithm generates some input variables $x_{ij}$ as arguments for the black-box function $f$. By utilizing the values of $f$, we may compute the values of the distribution ${{P}_{{{X}_{j}}}}(x_{ij})$, by Equation 6, and ${{F}_{{{X}_{j}}}}({{x}_{ij}})$ by Equation 7, for all $x_{ij}$.

In iteration $i$, the algorithm will have searched the space $\hat{A}\subset A$, within the limits $lb_{j},ub_{j}$. By definition, $P$ is the probability density (likelihood) that the optimum occurs in a region $\mathbf{x}\pm\mathbf{dx}$. Hence,

where $E[\cdot]$, denotes the expectation of a random variable or vector.

A variety of programming techniques were investigated, in order to implement the proposed method into a computer code. To keep the algorithm simple and reduce the computational time, we applied inverse transform sampling by keeping in each iteration $i$ the best function evaluations, and randomly sampling among them. This is equivalent to a kernel function that vanishes over the worst function evaluations, and distributing the probability mass to the best performing ones. Accordingly, similar programming techniques may be investigated in future works, within the suggested framework.

The Algorithm 2 represents a simple version of the supplementary code in the appendix, which may easily be programmed. The variable $\mathbf{opti\_evals}$ is a dynamic vector, containing all values of the objective function returned during the optimization history, until step $i$, and $\mathbf{x\_evals}$ is an $i\times n$ matrix, containing all the design vectors $\mathbf{x}_{1:i}$, corresponding to $\mathbf{opti\_evals}$. The integer $\alpha$ is a parameter regarding how many instances of the optimization history are kept in order to randomly sample among them; for example if $f_{e}=10^{4}$, we may select $\alpha=10^{2}$. In this variation of the code, the Inverse Transform Sampling is approximately implemented in line 7, by calculating the random variable $rr$, among the extrema of the vector $\mathbf{inp\_x}$, corresponding to the range where: for the $j^{th}$ dimension of all $\mathbf{x}_{1:i}$, the $\alpha$ best function values were returned by the black-box function $f$. The sought solution is the vector $\mathbf{opti\_x}$, and the mimimum attained value of the objective function $opti\_f$.

In this section, we present the results obtained by running the ITSO algorithm, as well as Adaptive Differential Evolution (rand 1 bin), Differential Evolution (rand 1 bin), and Differential Evolution (rand 2 bin) with and without radius limited, Compass Coordinate Direct Search, Probabilistic Descent Direct search, Random Search, Resampling Inheritance Memetic Search, Resampling Memetic Search, Separable Natural Evolution Strategies, Simultaneous Perturbation Stochastic Approximation, and Exponential Natural Evolution Strategies from the Julia Package BlackBoxOptim.jl [24], and Nelder-Mead, Particle Swarm, and Simulated Annealing from Optim.jl [25]. To demonstrate the performance of each optimizer in attaining the minimum, we firstly run the algorithm $r=10$ times, obtain $f_{k}(\mathbf{x}_{i})$ for $k=\{1,2,\dots,r\}$ and all iterations $i=\left\{1,2,\ldots,{{f}_{e}}\right\}$, and average the results

Then, we normalize the vector of obtained function evaluations $\mathbf{v}=\{\hat{f}(\mathbf{x}_{1}),\hat{f}(\mathbf{x}_{1}),\dots,\hat{f}(\mathbf{x}_{f_{e}})\}$ in the domain $\left[0,1\right]$ through

and finally use the optimization history $h$

where $m$ indicates the number of black-box functions used for the evaluation. Equation 22 was selected as a performance metric, in order to obtain a clear representation of the various optimizers utilized, as powers smaller than one (in our case $\frac{1}{10}$), have the property to magnify the attained values at the final steps of the optimization history. In Figure 3 the numerical experiments for $m=13$ functions are presented. Each line corresponds to the normalized average optimization history $h$ (Equation 22, for all functions, which were repeated 10 times). We may see a clear prevalence of the proposed framework, in terms of convergence performance, as expected by the theoretical investigation. The numerical experiments for the comparison with other optimizers, may be reproduced by running the file __run.jl.

In this work a novel approach was presented for the well known problem of finding the argument that minimizes a black-box, function or system. A vast volume of approximation algorithms have been proposed, mainly heuristic, such as genetic, evolutionary, particle swarm, as well as their variations. However, they stem from nature-inspired procedures, and hence their converge is investigated a-posteriori. Despite their efficiency, they are often deprecated by researchers, due to the lack of rigorous mathematical formulation, as well as complexity of implementation. To the contrary, the proposed algorithm, initiates its formulation from well established probabilistic definitions and theorems, and its implementation demands a few lines of computer code. Furthermore, the convergence properties were found stable, as a proof that the suggested framework attains the best possible solution in the fewest possible iterations. The numerical examples validate the theoretical results and may be reproduced by the provided computer code. We consider the suggested method as a powerful framework which may easily be adopted to the sought solution of any problem involving the minimization of a black-box function.