# MAC, a novel stochastic optimization method††thanks: This project is
supported by Hungarian National Research, Development, and Innovation Office
via grants FK134332 and K132109.

Attila László Nagy  Department of Stochastics, Institute of Mathematics,
Budapest University of Technology and Economics, Műegyetem rkp. 3., H-1111
Budapest, Hungary<EMAIL_ADDRESS> Goitom Simret Kidane  Institute of
Chemistry, ELTE Eötvös Loránd University, Budapest, Hungary  Tamás Turányi
Institute of Chemistry, ELTE Eötvös Loránd University, Budapest, Hungary
János Tóth  Institute of Chemistry, ELTE Eötvös Loránd University, Budapest,
Hungary  Department of Mathematical Analysis, Institute of Mathematics,
Budapest University of Technology and Economics, Műegyetem rkp. 3., H-1111
Budapest, Hungary<EMAIL_ADDRESS>
###### Abstract

A novel stochastic optimization method called MAC was suggested. The method is
based on the calculation of the objective function at several random points
and then an empirical expected value and an empirical covariance matrix are
calculated. The empirical expected value is proven to converge to the optimum
value of the problem. The MAC algorithm was encoded in Matlab and the code was
tested on 20 test problems. Its performance was compared with those of the
interior point method (Matlab name: fmincon), simplex, pattern search (PS),
simulated annealing (SA), particle swarm optimization (PSO), and genetic
algorithm (GA) methods. The MAC method failed two test functions and provided
inaccurate results on four other test functions. However, it provided accurate
results and required much less CPU time than the widely used optimization
methods on the other 14 test functions.

##  1 Introduction

Optimization is one of the central topics of both applied and theoretical
mathematics. It studies the problem of maximizing or minimizing outputs of
functions related to the selected parameters. Using efficient optimization,
one can in general increase efficiency and decrease losses. Many classes of
optimization problems have been studied in mathematics, computer science, and
other fields of science (e.g. see refs Hopfield and Tank, (1985); Krentel,
(1986); Gambella et al., (2021)). Over the years, a plethora of specialized
optimization areas emerged (e.g. convex optimization, stochastic programming,
combinatorial optimization, etc.) and various methods have been elaborated and
applied to specific problems. The algorithms can be classified into two groups
based on whether they use randomization (stochastic methods) or not
(deterministic methods). Stochastic optimization methods (e.g. Draper, (2002);
Spall, (2012); Zheng et al., (2015); Reddy et al., (2017)) search for an
optimal solution using the concept of probabilistic translation rules
(randomness). These algorithms are gaining popularity due to certain
properties which deterministic algorithms do not have. For example, stochastic
algorithms may suit to use for objective functions that have multiple local
optima whereby deterministic algorithms may get stuck.

An extremum (a maximum or minimum point) can be either global (the highest or
lowest function value within a region) or local (the highest or lowest value
in a neighborhood but not necessarily in the entire region). Finding a global
extremum, in general, proves to be an onerous task while trying to explore a
local extremum is usually less burdensome (see Press et al., (2007),
Zhigljavsky and Zilinskas, (2008), Weise, (2009) and Hendrix and Gazdag-Tóth,
(2010),). Several optimization methods have already been built into popular
program packages like Mathematica (Wolfram language) or Matlab.

The least-square function is a non-linear function of several variables
(called parameters), which is to be optimized on a convex and compact
parameter space. The least-square optimization can be regarded as one of the
most studied areas that are still investigated and applied in many different
areas of science and engineering (see Bard, (1974); Larrañaga and Lozano,
(2002); Seber and Wild, (2005)). The tasks of optimization are either
maximization or minimization, but they are trivially related to each other,
hence from now we are focusing on and formalizing only minimization problems.
We are interested in minimizing a real-valued continuous function defined over
a convex, compact domain. The solutions to these types of problems are often
based on finding the minimum of some kind of distance between the model
results and the corresponding measurements Bard, (1974); Gustavo et al.,
(2020).

In the applications of chemical kinetics to combustion systems, optimization
of rate parameters of reaction mechanisms is required for the efficient
interpretation of experimental data Goitom et al., (2022). The rate parameters
describe the temperature and pressure dependence of the rate coefficients of
the reaction steps, and their initial values are assigned based on direct
measurements of rate coefficients, theoretical calculations, or analogies to
similar reaction steps. The mechanism is then tested against indirect
experimental data, like measured laminar burning velocities and ignition delay
times, through an optimization process. This optimization is carried out, by
comparing the experimental and simulated data points using a nonlinear
objective function called an error function. This error function have
typically many local minima and the evaluation of the error function is
computationally expensive because it is based on the solution of large systems
of ODEs and PDEs. A search for efficient optimization methods applicable to
chemical kinetics optimizations is one of the motivations of this paper.

As our focus is on the optimization of a function over a domain, the two main
important issues to investigate are:

  * (i)

how many times we have to evaluate the error function and

  * (ii)

how many iterations are needed to arrive at or close to the extremum point?

In the present work, we describe a novel method that belongs to the family of
stochastic approximation methods. The method generates a converging series of
mean values and (empirical) covariance matrices, therefore it was named the
MAC method. It seems to have good properties from the above-enumerated
viewpoints.

Convergence of the method will be proved in Theorem 1. We have chosen a series
of benchmark functions found in Jamil and Yang, (2013); Mishra, (2006);
Rahnamayan et al., (2007); Shang and Qiu, (2006) and we have found that on
several test problems, the novel method requires much fewer function
evaluations, thus it is appropriate to solve problems when the function
evaluation is expensive. Also, it approaches close to the minimum in fewer
steps than several other methods.

The structure of the paper is as follows. Section 2.1 introduces the main
concepts and presents the algorithm. After that (Section 2.2) the mathematical
results are discussed. This is followed by the presentation of the numerical
results on benchmark functions in Section 3. We summarize the results and
outline the plans for the near future in Section 4.

##  2 Main concepts

This section is devoted to presenting our stochastic optimization algorithm
(MAC) in rigorous terms. After a preliminary discussion on the problem, we
introduce the main algorithm and later we discuss some results.

###  2.1 The algorithm

We are given a continuous, deterministic, real-valued function $Q$ of $d$
independent variables
($\mathbf{x}=(x_{1},x_{2},\ldots,x_{d})^{\top},d\in\mathbb{Z}^{+}$) defined on
a domain $\mathcal{D}\in\mathds{R}^{d}$. We assume that $\mathcal{D}$ is a
convex, compact (i.e. closed and bounded) set of $\mathds{R}^{d}$. Under these
terms, it follows that $Q$ is a bounded function and its (global and local)
minimum(s) and maximum(s) are attained over $\mathcal{D}$.

Our main goal is to find the extrema of $Q$ with a stochastic iterative
algorithm. Minimizing $Q$ over $\mathcal{D}$ is equivalent to maximizing $-Q$
above the same region. Hence, without loss of generality, we are targeting to
find the (global) minimum points of $Q$ on $\mathcal{D}$. Note that $Q$ is
equivalently called objective function throughout the article.

Under the above terms, let us also define the (non-empty) set $\mathcal{G}$ of
global minimum points of $Q$ over $\mathcal{D}$ as

$\mathcal{G}:=\left\\{\mathbf{x}\in\mathcal{D}:\,Q(\mathbf{x})\leq
Q(\mathbf{y})\,\mbox{ holds for all $\mathbf{y}\in\mathcal{D}$}\right\\}.$ (1)

Let us introduce the mass (or penalty) function $g$ that is associated with
the objective function $Q$ in such a way that for
$(\rho,\mathbf{x})\in\mathds{R}^{+}_{0}\times\mathds{R}^{d}$:

$g(\rho,\mathbf{x})=\exp\left(-\rho\,Q(\mathbf{x})\right).$ (2)

Note that $0<g(\rho,\mathbf{x})\leq 1$ holds for every
$\rho\in\mathds{R}^{+}_{0}$ and $\mathbf{x}\in\mathcal{D}$.

Next, we introduce the concept of ambient (apriori) distribution. This is
simply defined to be any fixed $\mathbf{\nu}$ absolutely continuous,
$d$-dimensional probability distribution such that:

  * •

the support of $\mathbf{\nu}$ is a convex, compact subset of $\mathds{R}^{d}$;

  * •

$\mathbf{\nu}$ is ’standard’, that is its expectation is $\mathbf{0}$ and its
covariance matrix is the identity matrix $\mathbf{I}$.

E.g. $\mathbf{\nu}$ can be chosen to be the standard $d$-dimensional normal
distribution (support of which is $\mathds{R}^{d}$) or it can be the uniform
distribution over the $d$-dimensional unit ’ball’ around the origin (that is
its support is $\\{\mathbf{x}:\left\|\mathbf{x}\right\|_{2}\leq 1\\}$).

Finally, let
$\alpha,\gamma:\mathbb{Z}^{+}_{0}\times\mathbb{Z}^{+}_{0}\rightarrow\mathbb{Z}^{+}_{0}$
be two non-decreasing, unbounded, double sequences of natural numbers, that is
for all $m\leq n,\,M\leq N\in\mathbb{Z}^{+}_{0}$: we have
$\alpha(m,M)\leq\alpha(n,N)$, $\gamma(m,M)\leq\gamma(n,N)$, and
$\sup_{l\in\mathbb{Z}^{+}}\alpha(l,N)=\sup_{L\in\mathbb{Z}^{+}}\alpha(n,L)=\sup_{l\in\mathbb{Z}^{+}}\gamma(l,N)=\sup_{L\in\mathbb{Z}^{+}}\gamma(n,L)=+\infty$
for all fixed $n,N\in\mathbb{Z}^{+}_{0}$. We also fix $\alpha(0,0)$=0.

Now, we are ready to present the pseudo-code of our stochastic algorithm.

  1. BEGIN

  2. 1.

let $\beta:\mathbb{Z}^{+}_{0}\rightarrow\mathbb{Z}^{+}_{0}$ be a fixed, non-
decreasing function of natural numbers such that $\beta(0)=0$ and let
$\delta>0$ be a real number; with a slight abuse of notation, we use the
shorthand notations below $\alpha_{n}:=\alpha(n,\beta(n))$ and
$\gamma_{n}:=\gamma(n,\beta(n))$ for $n\in\mathbb{Z}^{+}_{0}$.

  3. 2.

let $\mathbf{u}_{0}\in\mathcal{D}$, pick a positive definite matrix
$\mathbf{U}_{0}\in\mathds{R}^{d\times d}$, and fix an ambient (apriori)
distribution $\mathbf{\nu}$;

  4. 3.

$n:=1$;

    1. REPEAT

    2. i.

Take a completely independent sample
$\left(\xi_{i}\right)_{i=\alpha_{n-1}+1}^{\alpha_{n}}$ from the distribution
$\mathbf{\nu}$ and let:

$\mathbf{p}_{i}:=\mathbf{u}_{n-1}+\mathbf{U}_{n-1}\xi_{i},$

where $i=\alpha_{n-1}+1,\alpha_{n-1}+2,\ldots,\alpha_{n}$. (Note that
$\mathbf{p}_{i}$’s are already calculated in the previous step(s) for
$i=1,2,\ldots,\alpha_{n-1}$.) Finally, compute the weights:
$(\mathbf{p}_{i})_{i=\alpha_{n-1}}^{\alpha_{n}}$.

    3. ii.

Compute the empirical expected value and empirical covariance matrix with
those weights determined above, i.e.:

$\displaystyle\mathbf{u}_{n}$
$\displaystyle:=\sum_{i=1}^{\alpha_{n}}\frac{g(\gamma_{n},\mathbf{p}_{i})}{\sum_{j=1}^{\alpha_{n}}g(\gamma_{n},\mathbf{p}_{j})}\,\mathbf{p}_{i},$
(3) $\displaystyle\mathbf{U}_{n}$
$\displaystyle:=\left(\sum_{i=1}^{\alpha_{n}}\frac{g(\gamma_{n},\mathbf{p}_{i})}{\sum_{j=1}^{\alpha_{n}}g(\gamma_{n},\mathbf{p}_{j})}\,(\mathbf{p}_{i}-\mathbf{u}_{n})(\mathbf{p}_{i}-\mathbf{u}_{n})^{\top}\right)^{1/2},$
(4)

where recall the penalty function defined in Eq. 2.

    4. iii.

$n:=n+1$;

    5. UNTIL $\\{\left\|\mathbf{u}_{n}-\mathbf{u}_{n-1}\right\|_{2}<\delta\\}$.

  5. END

######  Remark.

The $\alpha_{n}$ and $\gamma_{n}$ parameters are also called the ’learning’
parameters. The more we tweak them, i.e. increase, the closer we might get to
the optimum (minimum) point.

Note that ’$\mathbf{U}_{n}^{2}$’ is a symmetric, positive semi-definite
matrix, hence its square root is well-defined and is real as it can be
decomposed as $\mathbf{O}\mathbf{\Lambda}\mathbf{O}^{\top}$ for some
orthogonal matrix $\mathbf{O}\in\mathds{R}^{d\times d}$, where
$\mathbf{\Lambda}$ is the diagonal matrix of the eigenvalues of
$\mathbf{U}_{n}^{2}$, hence
$\mathbf{U}_{n}=\mathbf{O}\mathbf{\Lambda}^{1/2}\mathbf{O}^{\top}$.

Mind that $\mathbf{u}_{n}\in\mathcal{D}$ holds for all $n\in\mathbb{Z}^{+}$
due to the convexity of the set $\mathcal{D}$, hence,
$\\{\mathbf{u}_{n}\\}_{n\in\mathbb{Z}^{+}}$ is a bounded sequence in
$\mathds{R}^{d}$. However, $\\{\mathbf{U}_{n}\\}_{n\in\mathbb{Z}^{+}}$ is only
bounded when the support of $\mathbf{\nu}$ is finite.

Now, the immediate question of whether the above method converges at all, more
precisely $\mathbf{u}_{n}$ converges to an extremum point, is going to be
answered affirmatively under some further conditions.

###  2.2 Results

The main theorem of this section concerns the proof of convergence for the
previously introduced algorithm under some further conditions on its
parameters.

######  Theorem 1.

Let $\mathbf{u}_{0}\in\mathcal{D}$ and $\mathbf{U}_{0}\in\mathds{R}^{d\times
d}$ be a column-vector and an arbitrarily chosen positive definite matrix,
respectively. Assume that the ambient (apriori) distribution $\nu$ is a.s.
bounded in such a way that for some $0<T<\frac{1}{4}$:

$\nu\left(\mathcal{S}\right)=1,\,\mbox{ where
$\,\mathcal{S}=\\{\mathbf{x}:\left\|\mathbf{x}\right\|_{2}^{2}\leq T\\}$}.$
(5)

Let
$\mathcal{A}:=\\{\mathbf{x}\in\mathcal{S}\,:\,\mathbf{u}_{0}+\mathbf{U}_{0}\mathbf{x}\in\mathcal{D}\\}$
and assume that $Q$ has a single minimum value over
$\mathbf{u}_{0}+\mathbf{U}_{0}\mathcal{A}$, i.e.:

$\displaystyle\operatorname*{arg\,min}_{\mathbf{x}\in\mathcal{A}}Q(\mathbf{u}_{0}+\mathbf{U}_{0}\mathbf{x})$
$\displaystyle=\\{\overline{\mathbf{x}}\\}.$

In other words,
$\mathcal{G}=\\{\mathbf{u}_{0}+\mathbf{U}_{0}\overline{\mathbf{x}}\\}$, where
recall 1. Furthermore, with a slight abuse of notation, let
$\alpha_{n}=N\cdot\tilde{\alpha}_{n}$ and
$\gamma_{n}=\gamma\cdot\tilde{\gamma}_{n}$, where $N\in\mathbb{Z}^{+}$,
$\gamma\in\mathds{R}^{+}$ are free parameters, and
$\tilde{\alpha},\tilde{\gamma}$ alone satisfy the requirements of the
algorithm listed above. Then, we have the following convergence results

$\displaystyle\lim_{(n,N,\gamma)\to+\infty}\mathbf{u}_{n}=\mathbf{u}_{0}+\mathbf{U}_{0}\overline{\mathbf{x}},$
(6)
$\displaystyle\lim_{(n,N,\gamma)\to+\infty}\mathbf{U}_{n}=\mathbf{0}\mathbf{0}^{\top}$
(7)

to hold such that for Eq. 6 we impose the condition
$\lim_{(N,\gamma)\rightarrow+\infty}N\cdot\exp\left(-\gamma\right)=0$, where
the above convergence is meant almost surely (in short: a.s.).

  * Proof

(Theorem 1). First, we prove the second convergence 7 (for empirical variance)
and then focus on the first one. First, let
$c=\mathchoice{{\hbox{$\displaystyle\sqrt{4T\,}$}\lower 0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.39828pt}}}{{\hbox{$\textstyle\sqrt{4T\,}$}\lower
0.4pt\hbox{\vrule
height=6.83331pt,depth=-5.39828pt}}}{{\hbox{$\scriptstyle\sqrt{4T\,}$}\lower
0.4pt\hbox{\vrule
height=4.78333pt,depth=-3.77881pt}}}{{\hbox{$\scriptscriptstyle\sqrt{4T\,}$}\lower
0.4pt\hbox{\vrule height=3.41666pt,depth=-2.69914pt}}}$ and note that by our
assumption 5: $0<c<1$ holds.

It is straightforward to see from Eq. 4 that

$\displaystyle\mathbf{U}_{n}^{2}=\mathbf{U}_{n-1}\left(\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}\mathbf{w}_{i}\mathbf{w}_{i}^{\top}\right)\mathbf{U}_{n-1},$

where we used the below short-hand notations:

$\displaystyle\mathbf{w}_{i}$
$\displaystyle=\xi_{i}-w_{i,\alpha_{n},\gamma_{n}}\xi_{i},$ $\displaystyle
w_{i,\alpha_{n},\gamma_{n}}$
$\displaystyle=\frac{g(\gamma_{n},\mathbf{p}_{i})}{\sum_{j=1}^{\alpha_{n}}g(\gamma_{n},\mathbf{p}_{j})}$

Recall that the random variables $\xi_{i}$’s are directly sampled from $\nu$,
hence $\xi_{i}\in\mathcal{S}$ holds for every $1\leq i\leq\alpha_{n}$. Then,
by 5, the following chain of inequality can easily be derived:

$\displaystyle\left\|\mathbf{w}_{i}\mathbf{w}_{i}^{\top}\right\|_{F}$
$\displaystyle\leq\sup_{\begin{subarray}{c}\kappa_{i}\in[0,1],\mathbf{x}_{i}\in\mathcal{S}\\\
\sum_{i=1}^{\alpha_{n}}\kappa_{i}=1,\left\|\mathbf{x}_{i}\right\|_{2}^{2}\leq
T\\\ (1\leq
i\leq\alpha_{n})\end{subarray}}\left\|\left(\mathbf{x}_{i}-\sum_{j=1}^{\alpha_{n}}\kappa_{i}\mathbf{x}_{i}\right)\left(\mathbf{x}_{i}^{\top}-\sum_{i=1}^{\alpha_{n}}\kappa_{i}\mathbf{x}_{i}^{\top}\right)\right\|_{F}$
$\displaystyle\leq\sup_{\mathbf{x}\in\mathcal{S}}\left\|2\mathbf{x}\right\|_{2}^{2}=4T=c^{2},$

holding for all $1\leq i\leq\alpha_{n}$, where $\left\|\cdot\right\|_{F}$
denoted the Frobenius norm. Now, from the above estimate it is clear that

$\displaystyle\left\|\mathbf{U}_{n}\right\|_{F}^{2}$
$\displaystyle=\operatorname{Tr}\mathbf{U}_{n}^{2}$
$\displaystyle=\operatorname{Tr}\mathbf{U}_{n-1}\left(\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}\mathbf{w}_{i}\mathbf{w}_{i}^{\top}\right)\mathbf{U}_{n-1}$
$\displaystyle=\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}\operatorname{Tr}\left(\mathbf{w}_{i}\mathbf{w}_{i}^{\top}\right)\mathbf{U}_{n-1}^{2}$
$\displaystyle\leq\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}\left|\left\langle\mathbf{w}_{i}\mathbf{w}_{i}^{\top},\mathbf{U}_{n-1}^{2}\right\rangle_{F}\right|$
$\displaystyle\leq\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}\left\|\mathbf{w}_{i}\mathbf{w}_{i}^{\top}\right\|_{F}\left\|\mathbf{U}_{n-1}^{2}\right\|_{F}$
$\displaystyle=\left\|\mathbf{U}_{n-1}\right\|_{F}^{2}\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}\left\|\mathbf{w}_{i}\right\|_{2}^{2}$
$\displaystyle\leq
c^{2}\left\|\mathbf{U}_{n-1}\right\|_{F}^{2}\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}(1-w_{i,\alpha_{n},\gamma_{n}})$
$\displaystyle=c^{2}\left\|\mathbf{U}_{n-1}\right\|_{F}^{2}\left(1-\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}^{2}\right)$
(8)

holds for every $n\in\mathbb{Z}^{+}$, where we used the cyclicity and
linearity of the trace, the Cauchy–Bunyakovsky–Schwarz inequality, and the
sub-multiplicative property of the Frobenius norm along the way.

Now, using the inequality just obtained in 8, an induction on $n$ leads to the
following estimates for $\mathbf{U}_{n}$:

$\displaystyle\left\|\mathbf{U}_{n}\right\|_{2}$
$\displaystyle\leq\left\|\mathbf{U}_{n}\right\|_{F}$ $\displaystyle\leq
c^{n-1}\left\|\mathbf{U}_{1}\right\|_{F}$ $\displaystyle\leq
c^{n}\left\|\mathbf{U}_{0}\right\|_{F}\mathchoice{{\hbox{$\displaystyle\sqrt{1-\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.09108pt}}}{{\hbox{$\textstyle\sqrt{1-\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=6.44444pt,depth=-5.09108pt}}}{{\hbox{$\scriptstyle\sqrt{1-\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}^{2}\,}$}\lower
0.4pt\hbox{\vrule
height=4.51111pt,depth=-3.56375pt}}}{{\hbox{$\scriptscriptstyle\sqrt{1-\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}^{2}\,}$}\lower
0.4pt\hbox{\vrule height=3.44165pt,depth=-2.71889pt}}}$ $\displaystyle\leq
c^{n}\left\|\mathbf{U}_{0}\right\|_{F}$ (9)

that hold for every $n\in\mathbb{Z}^{+}$, taking advantage of the subordinate
property of the Frobenius norm to the 2-norm and the fact that $\nu$ is being
finitely supported here. This then implies that
$\mathbf{U}_{n}\to\mathbf{0}\mathbf{0}^{\top}$ as $(n,N,\gamma)\to+\infty$
a.s. Since each norm is equivalent in finite dimension, $\mathbf{U}_{n}$
converges to $\mathbf{0}\mathbf{0}^{\top}$ a.s. in any (matrix) norm
completing the proof of Eq. 7.

Now, turning to the proof of Eq. 6, observe that using the above notation we
have

$\displaystyle\mathbf{u}_{n}=\mathbf{u}_{n-1}+\mathbf{U}_{n-1}\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}\xi_{i}$
(10)

to hold for every $n\in\mathbb{Z}^{+}$. As a consequence

$\displaystyle\left\|\mathbf{u}_{n}-\mathbf{u}_{n-1}\right\|_{2}$
$\displaystyle\leq\left\|\mathbf{U}_{n-1}\right\|_{2}\left\|\sum_{i=1}^{\alpha_{n}}w_{i,\alpha_{n},\gamma_{n}}\xi_{i}\right\|_{2}$
$\displaystyle\leq\left\|\mathbf{U}_{n-1}\right\|_{2}\sup_{\mathbf{x}\in\mathcal{S}}\left\|\mathbf{x}\right\|_{2}\leq
c\left\|\mathbf{U}_{n-1}\right\|_{F}\leq
c^{n}\left\|\mathbf{U}_{0}\right\|_{F},$

using the previously deduced inequality for $\mathbf{U}_{n}$ (Ineq. 9), and
the relation between the Frobenius- and the 2-norm of a matrix. With the above
inequality under the belt, it is not hard to see that $\mathbf{u}_{n}$ is
indeed a Cauchy sequence. Let $n,m\in\mathbb{Z}^{+}$ such that $n>m$, then it
holds that

$\displaystyle\left\|\mathbf{u}_{n}-\mathbf{u}_{m}\right\|_{2}$
$\displaystyle\leq\sum_{i=1}^{n-m}\left\|\mathbf{u}_{m+i}-\mathbf{u}_{m+i-1}\right\|_{2}$
$\displaystyle\leq\left\|\mathbf{U}_{0}\right\|_{F}\sum_{i=1}^{n-m}c^{m+i}\leq\left\|\mathbf{U}_{0}\right\|_{F}\frac{c^{m}}{1-c},$
(11)

which goes to zero as $(n,m)\rightarrow+\infty$. It then follows that there
exists a $\mathbf{u}_{\infty}$ random variable such that
$\mathbf{u}_{n}\rightarrow\mathbf{u}_{\infty}$ as $n\rightarrow+\infty$ a.s.,
and

$\displaystyle\left\|\mathbf{u}_{n}-\mathbf{u}_{\infty}\right\|_{2}\leq C\cdot
c^{n},$ (12)

where $C=\frac{1}{1-c}\left\|\mathbf{U}_{0}\right\|_{F}$ is an absolute
constant depending only on the choice of $\mathbf{U}_{0}$, that is no matter
how we have chosen $(N,\gamma)$.

Note that Eq. 10 can be further written as:

$\displaystyle\mathbf{u}_{n}=\mathbf{u}_{0}+\mathbf{U}_{0}\sum_{i=1}^{\alpha_{1}}w_{i,\alpha_{1},\gamma_{1}}\xi_{i}+\sum_{j=2}^{n}\mathbf{U}_{j-1}\sum_{i=1}^{\alpha_{j}}w_{i,\alpha_{j},\gamma_{j}}\xi_{i}\quad(n>0),$
(13)

and from this it is easy to see that for each fixed $n\in\mathbb{Z}^{+}$, the
sequence $(\mathbf{u}_{n})_{(N,\gamma)}$ as $(N,\gamma)\rightarrow+\infty$
converges a.s. using the strong law of large numbers. Hence, the limit
$\lim_{(n,N,\gamma)\rightarrow+\infty}\mathbf{u}_{n}$ exists a.s., taking
advantage of the uniform estimate pf Ineq. 12. Finally, using the fact that
$\sum_{i=1}^{\alpha_{1}}w_{i,\alpha_{1},\gamma_{1}}^{2}$ goes to $1$ as
$(N,\gamma)\rightarrow+\infty$ such that we have
$N\cdot\exp\left(-\gamma\right)\rightarrow 0$ in place by our assumption, and
by recalling Ineq. 9, we can conclude that the double sum (third term) of Eq.
13 converges to zero, whereas the rest approaches the global minimum.
$\blacksquare$

##  3 Benchmarking

The new method was tested on 18 artificial optimization test functions and the
results were compared with results obtained from a series of widely used
numerical optimization methods.

###  3.1 Setup for benchmarking

The newly proposed stochastic global optimization method, to be referred to as
the MAC method was implemented in the Matlab programming language. MAC refers
to empirical Mean And empirical Covariance matrix as it was based on a series
of calculations of empirical mean and empirical covariance matrix. Several
nonlinear uni-modal and multi-modal benchmark functions were selected to
compare the performance of the method with other well-known numerical
optimization methods. The functions given by an explicit formula that is used
for this benchmarking are outlined in Table 1, below. Recently, Layeb, (2022)
has found a few functions that are even harder to minimize that are included
in Table 1 (f9 to f18).

ID | Function Name | Dim | Domain
---|---|---|---
f1 | Ackley | 10 | $\left[-32.768\ 32.768\right]^{5}$
f2 | Cross-in-tray | 2 | $\left[-10\ 10\right]^{2}$
f3 | Rastrigin | 10 | $\left[-5.12\ 5.12\right]^{5}$
f4 | Rosenbrock | 10 | $\left[-5\ 10\right]^{5}$
f5 | RosenbrockSmall | 10 | $\left[-2.048\ 2.048\right]^{5}$
f6 | Rosenbrock scaled | 4 | $\left[0\ 1\right]^{4}$
f7 | Sphere | 20 | $\left[-5.12\ 5.12\right]^{5}$
f8 | Zakharov | 30 | $\left[-5\ 10\right]^{5}$
f9 | Layeb01 | 5 | $\left[-100\ 100\right]^{5}$
f10 | Layeb02 | 5 | $\left[-10\ 10\right]^{5}$
f11 | Layeb03 | 5 | $\left[-10\ 10\right]^{5}$
f12 | Layeb04 | 5 | $\left[-10\ 10\right]^{5}$
f13 | Layeb10 | 5 | $\left[-100\ 100\right]^{5}$
f14 | Layeb11 | 5 | $\left[-10\ 10\right]^{5}$
f15 | Layeb12 | 5 | $\left[-100\ 100\right]^{5}$
f16 | Layeb17 | 5 | $\left[-100\ 100\right]^{5}$
f17 | Layeb19 | 5 | $\left[-5\ 5\right]^{5}$
f18 | Layeb20 | 5 | $\left[-5\ 5\right]^{5}$
Table 1: Function ID, dimension, and domain of test functions

Details of the test functions (Table 1) are found in the following references
Jamil and Yang, (2013); Mishra, (2006); Rahnamayan et al., (2007); Shang and
Qiu, (2006) and Layeb, (2022). The MAC method is an iterative technique that
starts from some initial values of its hyper-parameters. These initial values
of the hyper-parameters depend on the user and problem type. There are two
important parameters that control the whole method (hence it has an effect on
its convergence, as well), they are:

  1. 1.

the sample size $\alpha_{n}$ that is generated at every step $n$, and

  2. 2.

another ’learning’ parameter denoted by $\gamma_{n}$.

As discussed in the pseudo-code, the MAC method has several tuning parameters.
In this implementation and testing of the MAC method, the following initial
values of these tuning parameters and their updating rules were used.

$\displaystyle N=4,\quad n=0$
$\displaystyle\gamma_{0}:=0.001,\quad\gamma_{n}:=2.8\times\gamma_{n-1}$
$\displaystyle\alpha_{0}:=1,\quad\alpha_{n}:=\alpha_{n-1}+n\times N.$

###  3.2 Benchmark results

The final calculated results of the different numerical methods and the MAC
method for the benchmark functions are given as tables in the following pages.
The interior point (Matlab name fmincon) and the simplex methods were used
from the MATLAB Optimization Toolbox, while the pattern search (PS), simulated
annealing (SA), particle swarm optimization (PSO), and genetic algorithm (GA)
methods were used from the MATLAB Global Optimization Toolbox. All the tested
methods were used with their default method parameters. In Table 2, the
minimum values are displayed, in Table 3, the number of function evaluations
required to obtain the minimum value, and in Table 4, the elapsed time (sec)
required are displayed. The calculations were performed on a Windows 10 PC,
Intel(R) Core(TM) i5-8250U CPU @ 1.80 GHz, RAM 8.00 GB.

Function | fmincon | simplex | PS | SA | PSO | GA | MAC | Exact
---|---|---|---|---|---|---|---|---
Ackley | $\mathrm{18.2}$ | 9.0 | $\mathrm{1.1E-4}$ | 15.6 | $\mathrm{2.2E-7}$ | $\mathrm{2.2E-5}$ | $\mathrm{8.0E-3}$ | 0
Cross-in-tray | -1.7906 | -1.7906 | -2.0626 | -2.0626 | -2.0626 | -2.0626 | -2.0626 | -2.0626
Rastrigin | fail | 74.6 | $\mathrm{5.9E-9}$ | 36.8 | 8.0 | $\mathrm{2.9E-7}$ | $\mathrm{2.7E-4}$ | 0
Rosenbrock | $\mathrm{1.7E-11}$ | $\mathrm{6.9E-9}$ | $\mathrm{2.6E-1}$ | $\mathrm{3.0E-1}$ | $\mathrm{8.3E-4}$ | $\mathrm{9.5E-2}$ | $\mathrm{8.5E-1}$ | 0
RosenbrockSmall | $\mathrm{6.4E-11}$ | $\mathrm{7.8E-24}$ | $\mathrm{1.7E-4}$ | 6.7 | 4.0 | $\mathrm{8.1E-4}$ | $\mathrm{1.6E-2}$ | 0
RosenbrockScaled | -1.0192 | -1.0192 | -1.0192 | -1.0191 | -1.0192 | -1.0191 | -1.0191 |
Sphere | $\mathrm{1.1E-23}$ | $\mathrm{2.5E-16}$ | $\mathrm{3.0E-11}$ | $\mathrm{4.3}$ | $\mathrm{1.3E-7}$ | $\mathrm{1.1E-5}$ | $\mathrm{2.3E-4}$ | 0
Zakharov | $\mathrm{5.4E-14}$ | fail | 74.6 | 9.5 | fail | $\mathrm{1.8E-4}$ | $\mathrm{3.7E-6}$ | 0
Layeb01 | 3.4 | fail | 0 | $\mathrm{3.8E-3}$ | fail | $\mathrm{6.0E-2}$ | 51.4 | 0
Layeb02 | $\mathrm{7.2E-1}$ | 0 | 0 | $\mathrm{6.0E-1}$ | $\mathrm{7.2E-1}$ | $\mathrm{5.8E-1}$ | 0 | 0
Layeb03 | $\mathrm{-8.7E-4}$ | $\mathrm{-2.3E-3}$ | $\mathrm{-1.3E-3}$ | $\mathrm{-3.8E-4}$ | fail | $\mathrm{-3.6E-3}$ | $\mathrm{-6.1E-4}$ | -4
Layeb04 | -31.631 | -11.561 | -31.631 | -7.761 | -31.631 | -31.631 | -23.6283 | -31.631
Layeb10 | fail | $\mathrm{1.3E-24}$ | 61.5 | 61.5 | fail | 75.3 | $\mathrm{2.9E-1}$ | 0
Layeb11 | fail | -2 | -4 | $\mathrm{-9.9E-1}$ | -2 | -3.0 | -3.6 | -4
Layeb12 | -12.3545 | -13.7137 | -14.8731 | -13.0683 | -13.2480 | -14.7952 | -14.785 | -14.873
Layeb17 | 4.0 | 16.8 | 0 | 57.6 | 4.0 | 4.0 | 4.3 | 0
Layeb19 | $\mathrm{2.5E-16}$ | $\mathrm{3.4E-30}$ | 0 | fail | $\mathrm{1.9E-13}$ | $\mathrm{2.2E-11}$ | $\mathrm{2.6E-6}$ | 0
Layeb20 | $\mathrm{1.4E-16}$ | $\mathrm{9.3E-28}$ | 0 | 27.7 | $\mathrm{2.0E-9}$ | $\mathrm{1.1E-6}$ | $\mathrm{2.8E-4}$ | 0
Table 2: The exact value versus the best value found by the different
optimization methods.

Function | fmincon | simplex | PS | SA | PSO | GA | MAC
---|---|---|---|---|---|---|---
Ackley | 163 | 1877 | 3096 | 5621 | 11100 | 5787 | 1723
Cross-in-tray | 24 | 191 | 234 | 1033 | 850 | 429 | 4810
Rastrigin | fail | 840 | 3406 | 4681 | 8700 | 11803 | 4615
Rosenbrock | 702 | 5652 | 10000 | 10001 | 230005 | 250563 | 4825
RosenbrockSmall | 755 | 4851 | 20000 | 8311 | 230050 | 470053 | 773
Rosenbrock Scaled | 400 | 567 | 1906 | 4511 | 2450 | 2920 | 2796
Sphere | 339 | 16323 | 11701 | 13500 | 10550 | 16362 | 5431
Zakharov | 2000 | fail | 20000 | 13441 | fail | 98988 | 2043
Layeb01 | 541 | fail | 321 | fail | fail | 32436 | 4831
Layeb02 | 78 | 60 | 76 | 8906 | 1300 | 9641 | 4810
Layeb03 | 252 | 504 | 722 | 3676 | fail | 9547 | 3233
Layeb04 | 644 | 555 | 704 | 3321 | 5300 | 6445 | 5431
Layeb10 | fail | 22856 | 191 | 2501 | fail | 7949 | 4810
Layeb11 | fail | 1880 | 153 | 5871 | 9050 | 5364 | 3233
Layeb12 | 200 | 979 | 266 | 3616 | 10600 | 9359 | 3713
Layeb17 | 480 | 2212 | 389 | 10306 | 3250 | 5834 | 4828
Layeb19 | 5005 | 1079 | 121 | fail | 2800 | 3578 | 2400
Layeb20 | 132 | 898 | 121 | 2626 | 3650 | 3390 | 4810
Table 3: The number of function evaluations until convergence. Function | fmincon | simplex | PS | SA | PSO | GA | MAC
---|---|---|---|---|---|---|---
Ackley | $\mathrm{7.4E-1}$ | $\mathrm{1.8E-1}$ | $\mathrm{2.4E-1}$ | $\mathrm{2.1}$ | $\mathrm{4.7E-1}$ | $\mathrm{1.6}$ | $\mathrm{2.7E-1}$
Cross-in-tray | $\mathrm{3.5E-2}$ | $\mathrm{7.1E-3}$ | $\mathrm{4.9E-2}$ | $\mathrm{3.3E-1}$ | $\mathrm{8.3E-2}$ | $\mathrm{3.1E-1}$ | $\mathrm{4.8E-1}$
Rastrigin | fail | $\mathrm{1.3E-1}$ | $\mathrm{5.2E-1}$ | $\mathrm{2.4}$ | $\mathrm{3.0E-1}$ | $\mathrm{1.0}$ | $\mathrm{3.0}$
Rosenbrock | $\mathrm{8.1E-1}$ | $\mathrm{3.4E-1}$ | $\mathrm{8.8E-1}$ | $\mathrm{11.2}$ | $\mathrm{3.3}$ | $\mathrm{13.6}$ | $\mathrm{9.2E-1}$
RosenbrockSmall | $\mathrm{8.0E-1}$ | $\mathrm{5.4E-1}$ | $\mathrm{2.2}$ | $\mathrm{4.1}$ | $\mathrm{4.1}$ | $\mathrm{16.4}$ | $\mathrm{1.2}$
RosenbrockScaled | $\mathrm{6.4E-1}$ | $\mathrm{7.5E-2}$ | $\mathrm{4.6}$ | $\mathrm{1.7}$ | $\mathrm{4.0E-1}$ | $\mathrm{7.1E-1}$ | $\mathrm{7.4E-1}$
Sphere | $\mathrm{5.6E-1}$ | $\mathrm{6.9E-1}$ | $\mathrm{6.0E-1}$ | $\mathrm{6.5}$ | $\mathrm{3.0E-1}$ | $\mathrm{1.4}$ | $\mathrm{1.5}$
Zakharov | $\mathrm{6.7E-1}$ | fail | $\mathrm{8.4E-1}$ | 8.0 | fail | 6.1 | $\mathrm{9.7E-2}$
Layeb01 | $\mathrm{7.0E-1}$ | fail | $\mathrm{2.6E-1}$ | fail | fail | $\mathrm{1.9}$ | $\mathrm{1.7E-1}$
Layeb02 | $\mathrm{4.5E-1}$ | $\mathrm{5.2E-2}$ | $\mathrm{3.3E-1}$ | 2.5 | $\mathrm{1.7E-1}$ | 1.1 | $\mathrm{3.9E-1}$
Layeb03 | $\mathrm{5.3E-1}$ | $\mathrm{4.2E-1}$ | $\mathrm{6.9E-1}$ | 1.2 | fail | $\mathrm{8.8E-1}$ | $\mathrm{2.3E-1}$
Layeb04 | $\mathrm{9.7E-1}$ | $\mathrm{6.2E-2}$ | $\mathrm{2.9E-1}$ | $\mathrm{9.9E-1}$ | $\mathrm{2.3E-1}$ | $\mathrm{6.9E-1}$ | $\mathrm{4.8E-1}$
Layeb10 | fail | 1.1 | $\mathrm{9.1E-1}$ | $\mathrm{4.0E-1}$ | fail | $\mathrm{4.1E-1}$ | $\mathrm{8.8E-2}$
Layeb11 | fail | $\mathrm{1.8E-1}$ | $\mathrm{6.8E-2}$ | 2.1 | $\mathrm{3.3E-1}$ | $\mathrm{7.6E-1}$ | $\mathrm{6.3E-2}$
Layeb12 | $\mathrm{6.3E-1}$ | $\mathrm{8.4E-2}$ | $\mathrm{2.6E-1}$ | 1.6 | $\mathrm{3.4E-1}$ | 1.0 | $\mathrm{1.3E-1}$
Layeb17 | $\mathrm{6.8E-1}$ | $\mathrm{1.4E-1}$ | $\mathrm{3.4E-1}$ | 2.9 | $\mathrm{7.4E-1}$ | $\mathrm{1.4E-1}$ | $\mathrm{1.2E-1}$
Layeb19 | 2.7 | $\mathrm{2.5E-1}$ | 1.5 | fail | 2.2 | 1.3 | $\mathrm{6.8E-1}$
Layeb20 | $\mathrm{6.0E-1}$ | $\mathrm{2.3E-1}$ | $\mathrm{3.3E-1}$ | $\mathrm{6.7E-1}$ | $\mathrm{6.9E-1}$ | $\mathrm{9.8E-1}$ | $\mathrm{9.4E-1}$
Table 4: Running time (sec) until convergence.

Ackley 10D function

Cross-in-tray 2D function

Rastrigin 10D function

Rosenbrock 10D function

Sphere 20D function

Zakharov 30D function

Figure 1: The objective function value corresponding to the number of function
evaluation for different benchmark functions using MAC and other numerical
optimization methods

From the above tabular results, we can say that, the MAC method in general has
good performance in most of the benchmark functions except in some of the most
complex test functions. In some of the most complex test functions (f9, f11
and f16), the MAC method failed to find an optimum value of the objective
function as all the other numerical methods were failed too. According to the
summaries of the above tables (final objective function value), the MAC method
outperformed almost all the other numerical optimization methods on the test
functions: Rastrigin (f3), Zakharov (f8), Layeb02 (f10), Layeb10 (f13) and
Layeb11 (f14). In addition to that, the MAC method outperformed the SA in all
the test functions except in the Layeb01 (f9) test function. In terms of the
number of function evaluations, in general, the MAC method required more
function evaluation than all the local optimization methods and less function
evaluation than the global optimization methods. The running time requirement
in general does not have a significant difference among all the methods and
all the test functions. However, in most of the test functions, the MAC method
ended up with less running time (sec) than the other global optimization
methods.

##  4 Summary and outlook

A ”good” algorithm (or heuristics) must be very ”efficient” in such a way as
to be applicable to real-life problems, such as e.g. model fitting, maximum
likelihood estimate, cost minimization, and performance maximization. What do
we mean by ”efficient”? On the one hand, we have to take into consideration
that $Q$ is probably expensive to evaluate, and we may or may not be able to
efficiently and/or accurately compute its gradient or the Hessian matrix. So
our task is rather to find a minimum point of $Q$ with as few evaluations as
possible. On the other hand, we must recognize what we are looking for, we
must have well-defined stopping criteria, convergence order, and a clear image
of the code.

New and new optimization methods are being published even nowadays. One reason
for it is that numerical optimization methods are widely used for solving a
series of optimization problems having arisen in science and technology. The
development of the MAC stochastic optimization method presented in this paper
was motivated by the task of the estimation of reaction rate parameters of
large chemical kinetics models. Such models consist of large systems of
ordinary or partial differential equations. The specialty of these
optimization problems is that the prior domain of uncertainty of the fitted
parameters is large and within this domain, the objective function contains a
large number of local minima. Also, the evaluation of the objective function
is very computationally expensive, since it includes the solution of many
large systems of differential equations. Such optimization tasks arise in many
fields including chemistry science and engineering, combustion science, and
systems biology (e.g. development of metabolic models). Our immediate plan is
to test the MAC method on a real-life chemical kinetics problem related to
combustion models.

##  5 Supplementary material

We are attaching as Supplementary material

  1. 1.

The Matlab code of the MAC method

  2. 2.

The benchmark functions.

## References

  * Bard, (1974) Bard, Y. (1974).  Nonlinear Parameter Estimation.  Academic Press, New York.
  * Draper, (2002) Draper, D. F. D. (2002).  Stochastic optimization: a review.  International Statistical Review, 70:315–349.
  * Gambella et al., (2021) Gambella, C., Ghaddar, B., and Naoum-Sawaya, J. (2021).  Optimization problems for machine learning: A survey.  European Journal of Operational Research, 290(3):807–828.
  * Goitom et al., (2022) Goitom, S. K., Papp, M., Kovács, M., Nagy, T., Zsély, I. G., Turányi, T., and Pál, L. (2022).  Efficient numerical methods for the optimisation of large kinetic reaction mechanisms.  Combustion Theory and Modelling, 26(6):1071–1097.
  * Gustavo et al., (2020) Gustavo, M. F., Székely, E., and Tóth, J. (2020).  Kinetic modeling of a consecutive enzyme-catalyzed enantioselective reaction in supercritical media.  ACS Omega, 5:26795–26806.
  * Hendrix and Gazdag-Tóth, (2010) Hendrix, E. M. T. and Gazdag-Tóth, B. (2010).  Introduction to Nonlinear and Global Optimization.  Springer, New York.
  * Hopfield and Tank, (1985) Hopfield, J. J. and Tank, D. W. (1985).  “neural” computation of decisions in optimization problems.  Biological cybernetics, 52(3):141–152.
  * Jamil and Yang, (2013) Jamil, M. and Yang, X.-S. (2013).  A literature survey of benchmark functions for global optimisation problems.  International Journal of Mathematical Modelling and Numerical Optimisation, 4(2):150–194.
  * Krentel, (1986) Krentel, M. W. (1986).  The complexity of optimization problems.  In Proceedings of the eighteenth annual ACM symposium on Theory of computing, pages 69–76.
  * Larrañaga and Lozano, (2002) Larrañaga, P. and Lozano, J. A. (2002).  Estimation of Distribution Algorithms.  Genetic Algorithms and Evolutionary Computation. Springer.
  * Layeb, (2022) Layeb, A. (2022).  New hard benchmark functions for global optimization.  arXiv preprint arXiv:2202.04606.
  * Mishra, (2006) Mishra, S. K. (2006).  Performance of repulsive particle swarm method in global optimization of some important test functions: A fortran program.  Available at SSRN 924339.
  * Press et al., (2007) Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (2007).  Numerical Recipes: The Art of Scientific Computing.  Cambridge University Press, Cambridge, 3rd edition.
  * Rahnamayan et al., (2007) Rahnamayan, S., Tizhoosh, H. R., and Salama, M. M. (2007).  A novel population initialization method for accelerating evolutionary algorithms.  Computers & Mathematics with Applications, 53(10):1605–1614.
  * Reddy et al., (2017) Reddy, S. S., Sandeep, V., and Jung, C.-M. (2017).  Review of stochastic optimization methods for smart grid.  Frontiers in Energy, 11:197–209.
  * Seber and Wild, (2005) Seber, G. A. F. and Wild, C. J. (2005).  Nonlinear Regression.  Wiley Interscience, Hoboken, New Jersey.
  * Shang and Qiu, (2006) Shang, Y.-W. and Qiu, Y.-H. (2006).  A note on the extended rosenbrock function.  Evolutionary Computation, 14(1):119–126.
  * Spall, (2012) Spall, J. C. (2012).  Stochastic optimization.  Handbook of computational statistics: Concepts and methods, pages 173–201.
  * Weise, (2009) Weise, T. (2009).  Global Optimization Algorithms—Theory and Application.  Self-Published Thomas Weise, 2nd edition.
  * Zheng et al., (2015) Zheng, Q. P., Wang, J., and Liu, A. L. (2015).  Stochastic optimization for unit commitment—a review.  IEEE Transactions on Power Systems, 30:1913–1924.
  * Zhigljavsky and Zilinskas, (2008) Zhigljavsky, A. and Zilinskas, A. (2008).  Stochastic Global Optimization, volume 9 of Springer Optimization and Its Applications.  Springer.