# Multi-Resolution Online Deterministic Annealing: A Hierarchical and
Progressive Learning Architecture

Christos N. Mavridis, , and John S. Baras The authors are with the Department
of Electrical and Computer Engineering and the Institute for Systems Research,
University of Maryland, College Park, USA. emails:{mavridis,
<EMAIL_ADDRESS>partially supported by the Defense Advanced Research
Projects Agency (DARPA) under Agreement No. HR00111990027, by ONR grant
N00014-17-1-2622, and by a grant from Northrop Grumman Corporation.

###### Abstract

Hierarchical learning algorithms that gradually approximate a solution to a
data-driven optimization problem are essential to decision-making systems,
especially under limitations on time and computational resources. In this
study, we introduce a general-purpose hierarchical learning architecture that
is based on the progressive partitioning of a possibly multi-resolution data
space. The optimal partition is gradually approximated by solving a sequence
of optimization sub-problems online, using gradient-free stochastic
approximation updates. As a consequence, a function approximation problem can
be defined within each subset of the partition and solved using the theory of
two-timescale stochastic approximation. This simulates an annealing process
and defines a robust and interpretable heuristic method to gradually increase
the complexity of the learning architecture in a task-agnostic manner, giving
emphasis to regions of the data space that are considered more important
according to a predefined criterion. Finally, by imposing a tree structure in
the progression of the partitions, we provide a means to incorporate potential
multi-resolution structure of the data space into this approach, significantly
reducing its complexity, while introducing hierarchical variable-rate feature
extraction properties similar to certain classes of deep learning
architectures. Asymptotic convergence analysis and experimental results are
provided for supervised and unsupervised learning problems.

###### Index Terms:

Hierarchical Learning, Progressive Learning, Online Deterministic Annealing,
Multi-resolution Learning

##  I Introduction

Learning from observations is pivotal to autonomous decision-making and
communication systems. Mathematically, such learning problems are often
formulated as constrained stochastic optimization problems: given realizations
of a random variable $X\in S$ representing the observations, an optimal
parameter vector $\theta\in\Theta$ is to be found such that a well-defined
error measure between an unknown function $f(X)\in\mathcal{F}$ and a learning
model $\hat{f}(X,\theta)\in\mathcal{F}$, parameterized by $\theta$, is
minimized under potentially additional constraints. However, the solution of
such problems over the entire domain $S$ often requires the learning model
$\hat{f}(X,\theta)$ to be particularly complex, making the estimation of
$\theta$ costly, and raising issues with respect to phenomena such as over-
fitting, generalization, and robustness, connected by an underlying trade-off
between complexity and performance [1]. As a result, the ability to gradually
approximate a solution to these problems is essential to decision-making
systems that often operate in real-time and under limitations in memory and
computational resources.

Current deep learning methods have made progress towards the construction of a
hierarchical representation of the data space [2, 3, 4, 5]. However, such
approaches do not necessarily satisfy the above description of hierarchical
learning, since they typically use overly complex models over the entire data
space $S$, which comes in the expense of time, energy, data, memory, and
computational resources [6, 7]. In this work, we are mainly focusing on a
framework for hierarchical progressive learning and data representation, where
a gradually growing and hierarchically structured set of learning models is
used for function approximation. We consider a prototype-based learning
framework where, given random observations of $X\in S$, a set of prototypes
$\left\\{\mu_{i}\right\\}\in S$ (also called codevectors or neurons) are
scattered in the data space $S$ to encode subsets/regions
$\left\\{S_{i}\right\\}$ that form a partition of $S$ [8]. This adheres to the
principles of vector quantization for signal compression [9]. In this regard,
a knowledge representation can be defined as the set of codevectors
$\left\\{\mu_{i}\in S\right\\}$ that induce a structured partition
$\left\\{S_{i}\right\\}$ of the data space $S$, along with a set of local
learning models $\hat{f}(x,\theta_{i})$ associated with each region $S_{i}$,
parameterized by their own set of parameters $\theta_{i}$. A structured
representation like this allows, among other things, to locate specific
regions of the space that the algorithm needs to approximate in greater
detail, according to the problem at hand and the designer’s requirements. This
results in adaptively allocating more resources only in the subsets of the
data space that are needed, and provides benefits in terms of time, memory,
and model complexity. Moreover, learning with local models that take advantage
of the differences in the underlying distribution of the data space provides a
means to understand certain properties of the data space itself, i.e., this is
an interpretable learning approach [10]. An illustration of this framework is
given in Fig. 1.

(a) Classical regression problem. asdf asdf asdf asdf

(b) Combined problem of partitioning and function approximation.

(c) Tree-structured partitioning and function approximation.

Figure 1: Comparison of the classical regression problem over the entire
domain $S$ with the problem of combined partitioning and regression within
each subset of the partition. Here the input $x\in S$ is a random variable and
the function $f(x)$ is to be estimated over $S$ by (a) a single learning model
$\hat{f}(x,\theta)$, and (b)-(c) a set of
$\left\\{\hat{f}(x,\theta_{i})\right\\}$ defined in each region $S_{i}$, where
$\left\\{S_{i}\right\\}$ is a partition of S to be estimated as well.

Regarding the learning process, we are interested in algorithms that are able
to simultaneously solve both the problems of partitioning and function
approximation, given online (e.g., real-time) observations. This is of great
importance in many applications, and especially in the scope of learning
algorithms for inference and control in general cyber-physical systems [11,
12, 13]. To construct a sequence of partitions with increasing number of
subsets we build upon the notion of Online Deterministic Annealing [14] and
define a series of soft-clustering optimization problems:

$\displaystyle\min_{\left\\{\mu_{i}\right\\}}~{}F_{\lambda}(X,Q):=(1-\lambda)D(X,Q)-\lambda
H(X,Q),$

parameterized by a Lagrange coefficient $\lambda\in[0,1]$ controlling the
trade-off between minimizing an average distortion measure
$D(X,Q):=\mathbb{E}\left[d(X,Q)\right]$, for an appropriately defined
dissimilarity measure $d$, and maximizing the Shannon entropy $H(X,Q)$, with
$H(X,Q):=\mathbb{E}\left[-\log p(X,Q)\right]$. The novelty of the approach
lies in the introduction of $Q$ as a random variable described by the
association probabilities $p(\mu_{i}|X=x)$ that represents the probability of
a data point $x$ to belong to the subset $S_{i}:=\left\\{x\in
S:i=\operatorname*{arg\,min}_{j}d(x,\mu_{j})\right\\}$. Once the joint
probability space of $(X,Q)$ is defined, successively solving the optimization
problems $\min_{\left\\{\mu_{i}\right\\}}~{}F_{\lambda}(X,Q)$ for decreasing
values of $\lambda$, leads in a series of bifurcation phenomena when the
cardinality of the set of codevectors $\left\\{\mu_{i}\right\\}$ increases,
resembling an annealing process that introduces inherent robustness and
regularization properties [14, 15].

An important property of this approach, initially shown in [14], is that the
optimization problems $\min_{\left\\{\mu_{i}\right\\}}~{}F_{\lambda}(X,Q)$ can
be solved online, using gradient-free stochastic approximation updates [16],
as long as the measure $d$ belongs to the family of Bregman divergences,
information-theoretic dissimilarity measures that include, among others, the
widely used squared Euclidean distance and Kullback-Leibler divergence [17,
18]. We exploit the fact that a stochastic approximation algorithm can be used
as a training rule for constructing the partition $\left\\{S_{i}\right\\}$, to
build a framework that simultaneously trains the learning models
$\left\\{\hat{f}(x,\theta_{i})\right\\}$ defined in each region $S_{i}$. In
particular, according to the theory of two-timescale stochastic approximation
[16], we define two stochastic approximation algorithms that run at the same
time and with the same observations but with different stepsize schedules that
define a fast and a slow learing process. In our case the slow process
approximates the parameters $\left\\{\mu_{i}\right\\}$ and as a result the
partition $\left\\{S_{i}\right\\}$, and the fast process executes a function
approximation algorithm within each $S_{i}$ to find the optimal parameters
$\theta_{i}$ for the learning model $\hat{f}(x,\theta_{i})$.

Finally, we further extend this approach by incorporating structural
constraints in the construction of the partition $\left\\{S_{i}\right\\}$. In
particular, by imposing a non-binary tree structure in the growing set of the
parameters $\left\\{\mu_{i}\right\\}$, we show that we can both (a) greatly
reduce the quadratic (in the number of parameters $\mu_{i}$) complexity of the
approach, and (b) construct a hierarchical and progressively growing tree-
structured partition where each layer of the tree is trained using different
resolution representation of the data space, according to an independent
multi-resolution analysis. While this is a general framework for multi-
resolution learning, we show that, in the case when convolution-based multi-
resolution features are used, the proposed architecture shares similarities
with deep learning approaches such as Deep Convolutional Networks [2] and
Scattering Convolutional Networks [19]. Lastly, we provide asymptotic
convergence analysis of the proposed learning architecture and experimental
results to illustrate its properties in clustering, classification, and
regression applications.

The paper is organized as follows: Section II introduces the Online
Deterministic Annealing framework for progressive partitioning along with a
mathematical analysis of its properties. Section III develops the two-
timescale framework for combined partitioning and function approximation.
Section IV handles the problem of classification in two different approaches.
Section V extends the general model by incorporating tree-structure
constraints and multi-resolution representation of the data space. Finally,
Section VI illustrates experimental results, and Section VII concludes the
paper.

##  II Online Deterministic Annealing for Progressive Partitioning

In this section we provide a comprehensive review of the online deterministic
annealing approach introduced in [14] and [11], as well as additional
analytical results and insights that will be used in Sections III, IV, and V,
to construct the proposed hierarchical learning architecture.

We start our analysis with the case of unsupervised learning, where
partitioning a space $S$ is equivalent to the problem of clustering and
density estimation. In this context, the observations (data) are independent
realization of a random variable $X:\Omega\rightarrow S$ defined in a
probability space $\left(\Omega,\mathcal{F},\mathbb{P}\right)$, where
$S\subseteq\mathbb{R}^{d}$ is the observation space (data space). In a
prototype-based learning approach one defines a similarity measure
$d:S\rightarrow ri(S)$, and a set of $K$ prototypes/codevectors
$\mu:=\left\\{\mu_{i}\right\\}_{i=1}^{K}$, $\mu_{i}\in ri(S)$, that define a
partition $\left\\{S_{i}:x\in
S:i=\operatorname*{arg\,min}_{j}d(x,\mu_{j})\right\\}$ such that the following
average distortion measure is minimized:

$\min_{\mu}~{}J(\mu):=\mathbb{E}\left[\sum_{i}\mathds{1}_{\left[X\in
S_{i}\right]}d(X,\mu_{i})\right]$ (1)

Here $ri(S)$ represents the relative interior of $S$, and $\mathds{1}_{A}$ is
the indication function of an event $A$. The similarity measure as well as the
number of prototypes $K$ are predefined designer parameters. This process is
equivalent to finding the most suitable model out of a set of $K$ local
constant models, and results in a piecewise-constant approximation of the data
space $S$. This representation has been used for clustering in vector
quantization applications [20, 9], and, in the limit $K\rightarrow\infty$, can
be used for density estimation.

To construct a method that progressively increases the number of prototypes
$K$, we adopt a probabilistic approach similar to [11, 14], and define a
discrete random variable $Q:S\rightarrow ri(S)$ such that (1) takes the form

$\displaystyle\min_{\mu}~{}D(\mu)$
$\displaystyle:=\mathbb{E}\left[d\left(X,Q\right)\right]$ (2)
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[d(X,Q)|X\right]\right]$
$\displaystyle=\int p(x)\sum_{i}p(\mu_{i}|x)d(x,\mu_{i})~{}dx$

Notice that $Q$ is completely described by the association probabilities
$\\{p(\mu_{i}|x):=\mathbb{P}[Q=\mu_{i}|X=x]\\}$, $\forall i$. This is now a
problem of finding both the locations $\left\\{\mu_{i}\right\\}$ and the
association probabilities $\left\\{p(\mu_{i}|x)\right\\}$. Therefore this is a
more general problem than (1), where it is subtly assumed that
$p(\mu_{i}|x)=\mathds{1}_{\left[x\in S_{i}\right]}$.

The definition of the random variable $Q$ allows us to constraint the
distribution of $(X,Q)$ by maximizing the entropy:

$\displaystyle H(\mu)$ $\displaystyle:=\mathbb{E}\left[-\log
P(X,Q)\right]=H(X)+H(Q|X)$ (3) $\displaystyle=H(X)-\int
p(x)\sum_{i}p(\mu_{i}|x)\log p(\mu_{i}|x)~{}dx,$

at different levels. This is essentially a realization of the Jaynes’s maximum
entropy principle [21]. We formulate this multi-objective optimization as the
minimization of the Lagrangian

$\min_{\mu}F_{\lambda}(\mu):=(1-\lambda)D(\mu)-\lambda H(\mu)$ (4)

where $\lambda\in[0,1)$ acts as a Lagrange multiplier. The term
$T:=\frac{\lambda}{1-\lambda},\ \lambda\in[0,1)$ can be seen as a temperature
coefficient in a deterministic annealing process [14]. In this regard, this
approach follows from the Online Deterministic Annealing (ODA) algorithm in
[14], and its offline predecessor [15]. Equation (4) represents the
scalarization method for trade-off analysis between two performance metrics,
one related to performance, and one to generalization. The entropy $H$, acts
as a regularization term, and is given progressively less weight as $\lambda$
(resp. $T$) decreases. For large values of $\lambda\rightarrow 1$ (resp.
$T\rightarrow\infty$) we essentially maximize the entropy, and as $\lambda$
(resp. $T$) is lowered, we transition from one Pareto point to another in a
naturally occurring direction that resembles an annealing process.

In the remaining section, we will (i) derive an analytical solution of the
optimization problem (4) and a recursive gradient-free training rule to
approximate it online, (ii) show that the number of unique locations
$\left\\{\mu_{i}\right\\}$ is finite for $\lambda>0$ and increases as
$\lambda$ decreases beyond certain critical values with respect to a
bifurcation phenomenon, and (iii) analyze the asymptotic behavior and
complexity of this approach.

###  II-A Solving the Optimization Problem

As in the case of standard vector quantization algorithms, we will minimize
$F_{\lambda}$ in (4) by successively minimizing it first respect to the
association probabilities $\left\\{p(\mu_{i}|x)\right\\}$, and then with
respect to the codevector locations $\mu$. The following lemma provides the
solution of minimizing $F_{\lambda}$ with respect to the association
probabilities $p(\mu_{i}|x)$:

######  Lemma 1.

The solution of the optimization problem

$\displaystyle F_{\lambda}^{*}(\mu)$
$\displaystyle:=\min_{\left\\{p(\mu_{i}|x)\right\\}}F_{\lambda}(\mu)$ (5) s.t.
$\displaystyle\sum_{i}p(\mu_{i}|x)=1$

is given by the Gibbs distributions

$p^{*}(\mu_{i}|x)=\frac{e^{-\frac{1-\lambda}{\lambda}d(x,\mu_{i})}}{\sum_{j}e^{-\frac{1-\lambda}{\lambda}d(x,\mu_{j})}},~{}\forall
x\in S$ (6)

###### Proof.

See Appendix A. ∎

In order to minimize $F_{\lambda}^{*}(\mu)$ with respect to the codevector
locations $\mu$ we observe that

$\displaystyle\frac{d}{d\mu}F_{\lambda}^{*}(\mu)=\int
p(x)\sum_{i}(1-\lambda)\frac{d}{d\mu}\left(p^{*}(\mu_{i}|x)d_{\phi}(x,\mu_{i})\right)$
$\displaystyle\quad\quad\quad+\lambda\frac{d}{d\mu}\left(p^{*}(\mu_{i}|x)\log
p^{*}(\mu_{i}|x)\right)~{}dx$
$\displaystyle\phantom{\frac{d}{d\mu}F_{\lambda}^{*}(\mu)}=\int
p(x)\sum_{i}(1-\lambda)\frac{d}{d\mu}p^{*}(\mu_{i}|x)d_{\phi}(x,\mu_{i})$
$\displaystyle\quad\quad\quad+(1-\lambda)p^{*}(\mu_{i}|x)\frac{d}{d\mu}d_{\phi}(x,\mu_{i})+\lambda\frac{d}{d\mu}p^{*}(\mu_{i}|x)$
$\displaystyle\quad\quad\quad+\lambda\frac{d}{d\mu}p^{*}(\mu_{i}|x)\log
p^{*}(\mu_{i}|x)~{}dx$
$\displaystyle\phantom{\frac{d}{d\mu}F_{\lambda}^{*}(\mu)}=\int
p(x)\sum_{i}(1-\lambda)p^{*}(\mu_{i}|x)\frac{d}{d\mu}d_{\phi}(x,\mu_{i})$
$\displaystyle+\lambda\frac{d}{d\mu}p^{*}(\mu_{i}|x)-\lambda\frac{d}{d\mu}p^{*}(\mu_{i}|x)\sum_{j}e^{-\frac{1-\lambda}{\lambda}d_{\phi}(x,\mu_{j})}~{}dx$
$\displaystyle\phantom{\frac{d}{d\mu}F_{\lambda}^{*}(\mu)}=\sum_{i}\int
p(x)p^{*}(\mu_{i}|x)\frac{d}{d\mu_{i}}d(x,\mu_{i})~{}dx$

such that

$\frac{d}{d\mu}F_{\lambda}^{*}(\mu)=0\implies\sum_{i}\int
p(x)p^{*}(\mu_{i}|x)\frac{d}{d\mu_{i}}d(x,\mu_{i})~{}dx=0$ (7)

where we have used (6), direct differentiation, and
$\sum_{i}\frac{d}{d\mu}p^{*}(\mu_{i}|x)=\frac{d}{d\mu}\sum_{i}p^{*}(\mu_{i}|x)=0$.
In the following section, we show that (7) has an easy to compute closed form
solution if the dissimilarity measure $d$ belongs to the family of Bregman
divergences.

###  II-B Bregman Divergences as Dissimilarity Measures

The proximity measure $d$ can be generalized to dissimilarity measures
inspired by information theory and statistical analysis. In particular, the
family of Bregman divergences can offer numerous advantages in learning
applications compared to the Euclidean distance alone [17].

######  Definition 1 (Bregman Divergence).

Let $\phi:S\rightarrow\mathbb{R}$, be a strictly convex function defined on a
vector space $S\subseteq\mathbb{R}^{d}$ such that $\phi$ is twice
F-differentiable on $S$. The Bregman divergence $d_{\phi}:H\times
S\rightarrow\left[0,\infty\right)$ is defined as:

$\displaystyle
d_{\phi}\left(x,\mu\right)=\phi\left(x\right)-\phi\left(\mu\right)-\frac{\partial\phi}{\partial\mu}\left(\mu\right)\left(x-\mu\right),$

where $x,\mu\in S$, and the continuous linear map
$\frac{\partial\phi}{\partial\mu}\left(\mu\right):S\rightarrow\mathbb{R}$ is
the Fréchet derivative of $\phi$ at $\mu$.

The derivative of $d_{\phi}$ with respect to the second argument can be
written as

$\displaystyle\frac{\partial
d_{\phi}}{\partial\mu}(x,\mu)=-\frac{\partial^{2}\phi(\mu)}{\partial\mu^{2}}(x-\mu)=-\left<\nabla^{2}\phi(\mu),(x-\mu)\right>$
(8)

which leads to the following theorem showing that if $d$ is a Bregman
divergence, the solution to the second optimization step (7) can be
analytically computed in a convenient centroid form:

######  Theorem 2.

A sufficient condition for the solution of the optimization problem

$\min_{\mu}F_{\lambda}^{*}(\mu)$ (9)

where $F_{\lambda}^{*}(\mu)$ is defined in (5), is given by

$\mu_{i}^{*}=\mathbb{E}\left[X|\mu_{i}\right]=\frac{\int
xp(x)p^{*}(\mu_{i}|x)~{}dx}{p^{*}(\mu_{i})}$ (10)

if $d:=d_{\phi}$ is a Bregman divergence for some function $\phi$ that
satisfies Definition 1.

###### Proof.

Given (8), (7) becomes

$\int(x-\mu_{i})p(x)p^{*}(\mu_{i}|x)~{}dx=0$ (11)

which is equivalent to (10) since $\int
p(x)p^{*}(\mu_{i}|x)~{}dx=p^{*}(\mu_{i})$. ∎

As a final note, the family of Bregman divergences includes two notable
examples. The first is the widely used squared Euclidean distance
$d_{\phi}(x,\mu)=\|x-\mu\|^{2}$ ($\phi(x)=\left<x,x\right>,\
x\in\mathbb{R}^{d}$), and the second is the generalized Kullback-Leibler
divergence $d_{\phi}(x,\mu)=\left<x,\log
x-\log\mu\right>-\left<\mathds{1},x-\mu\right>$ ($\phi(x)=\left<x,\log
x\right>,\ x\in\mathbb{R}_{++}^{d}$).

###  II-C The Online Learning Rule

In an offline approach, the approximation of the conditional expectation
$\mathbb{E}\left[X|\mu_{i}\right]$ is computed by the sample mean of the data
points weighted by their association probabilities $p(\mu_{i}|x)$ [15]. To
define an online training rule for the deterministic annealing framework, a
stochastic approximation algorithm can be formulated [14] to recursively
estimate $\mathbb{E}\left[X|\mu_{i}\right]$ directly. The following theorem
follows directly from [11] and provides an online learning rule that solves
the optimization problem of (9).

######  Theorem 3.

Let $\left\\{x_{n}\right\\}$ be a sequence of independent realizations of $X$.
Then $\mu_{i}(n)$, defined by the online training rule

$\begin{cases}\rho_{i}(n+1)&=\rho_{i}(n)+\alpha(n)\left[\hat{p}(\mu_{i}|x_{n})-\rho_{i}(n)\right]\\\
\sigma_{i}(n+1)&=\sigma_{i}(n)+\alpha(n)\left[x_{n}\hat{p}(\mu_{i}|x_{n})-\sigma_{i}(n)\right]\end{cases}$
(12)

where $\sum_{n}\alpha(n)=\infty$, $\sum_{n}\alpha^{2}(n)<\infty$, and the
quantities $\hat{p}(\mu_{i}|x_{n})$ and $\mu_{i}(n)$ are recursively updated
as follows:

$\displaystyle\mu_{i}(n)=\frac{\sigma_{i}(n)}{\rho_{i}(n)},\quad\hat{p}(\mu_{i}|x_{n})=\frac{\rho_{i}(n)e^{-\frac{1-\lambda}{\lambda}d(x_{n},\mu_{i}(n))}}{\sum_{i}\rho_{i}(n)e^{-\frac{1-\lambda}{\lambda}d(x_{n},\mu_{i}(n))}}$
(13)

converges almost surely to a locally asymptotically stable solution of the
optimization (9), as $n\rightarrow\infty$.

###### Proof.

See Appendix B. ∎

The learning rule (12), (13) is a stochastic approximation algorithm [16]. In
the limit $\lambda\rightarrow 0$, it results in a consistent density estimator
according to the following theorem:

######  Theorem 4.

In the limit $\lambda\rightarrow 0$, and as the number of observed samples
$\left\\{x_{n}\right\\}$ goes to infinity, i.e., $n\rightarrow\infty$, the
learning algorithm based on (12), (13), results in a codebook $\mu$ that
constructs a consistent density estimator with
$\hat{p}(x)=\frac{\sum_{i}\mathds{1}_{\left[x\in S_{i}\right]}}{nVol(S_{i})}$,
where $S_{i}=\left\\{x\in
S:i=\operatorname*{arg\,min}\limits_{j}~{}d(x,\mu_{j})\right\\}$.

###### Proof.

See Appendix C. ∎

This means that as $\lambda\rightarrow 0$, the representation of the random
variable $X\in S$ by the codevectors $\mu$ becomes all the more accurate in
$S$, according to the underlying probability density $p(x)$. Regarding
clustering, the nearest-neighbor rule can be used to partition the space $S$
in Voronoi cells $S_{i}=\left\\{x\in
S:i=\operatorname*{arg\,min}_{j}~{}d_{\phi}(x,\mu_{j})\right\\}$.

######  Remark 1.

Notice that we can express the dynamics of the codevector parameters
$\mu_{i}(n)$ directly as:

$\displaystyle\mu_{i}(n+1)$
$\displaystyle=\frac{\alpha(n)}{\rho_{i}(n)}\bigg{[}\frac{\sigma_{i}(n+1)}{\rho_{i}(n+1)}(\rho_{i}(n)-\hat{p}(\mu_{i}|x_{n}))$
(14)
$\displaystyle\quad\quad\quad\quad\quad+(x_{n}\hat{p}(\mu_{i}|x_{n})-\sigma_{i}(n))\bigg{]}$

where the recursive updates take place for every codevector $\mu_{i}$
sequentially. This is a discrete-time dynamical system that presents
bifurcation phenomena with respect to the parameter $\lambda$, i.e., the
number of equilibria of this system changes with respect to the value
$\lambda$ which is hidden inside the term $\hat{p}(\mu_{i}|x_{n})$ in (13).
According to this phenomenon, the number of distinct values of $\mu_{i}$ is
finite, and the updates need only be taken with respect to these values that
we call “effective codevectors”. This is discussed in Section II-D.

###  II-D Bifurcation Phenomena

So far, we have assumed a countably infinite set of codevectors. In this
section we will show that the unique values of the set
$\left\\{\mu_{i}\right\\}$ that solves (4), form a finite set $K(\lambda)$ of
values that we will refer to as “effective codevectors” throughout this paper.
In other words, both the number and the locations of the codevectors depend on
the value of $\lambda$ (resp. the value of the temperature parameter $T$).
These effective codevectors are the only values that an algorithmic
implementation will need to store in memory and update.

First, notice that when $\lambda\rightarrow 1$ (resp. $T\rightarrow\infty$)
equation (6) yields uniform association probabilities
$p(\mu_{i}|x)=p(\mu_{j}|x),\ \forall i,j,\forall x$. As a result of (7), all
codevectors are located at the same point:

$\displaystyle\mu_{i}=\mathbb{E}\left[X\right],\ \forall i$

which means that there is one unique effective codevector given by
$\mathbb{E}\left[X\right]$.

As $\lambda$ is lowered below a critical value, a bifurcation phenomenon
occurs, when the number of effective codevectors increases. Mathematically,
this occurs when the existing solution $\mu^{*}$ given by (10) is no longer
the minimum of the free energy $F^{*}$, as $\lambda$ (resp. the temperature
$T$) crosses a critical value. Following principles from variational calculus,
we can rewrite the necessary condition for optimality (7) as

$\frac{d}{d\epsilon}F^{*}(\mu+\epsilon\psi)|_{\epsilon=0}=0$ (15)

with the second order condition being

$\frac{d^{2}}{d\epsilon^{2}}F^{*}(\left\\{\mu+\epsilon\psi\right\\})|_{\epsilon=0}\geq
0$ (16)

for all choices of finite perturbations $\left\\{\psi\right\\}$. Here we
denote by $\left\\{y:=\mu+\epsilon\psi\right\\}$ a perturbed codebook, where
$\psi$ are perturbation vectors applied to the codevectors $\mu$, and
$\epsilon\geq 0$ is used to scale the magnitude of the perturbation.
Bifurcation occurs when equality is achieved in (16) and hence the minimum is
no longer stable111For simplicity we ignore higher order derivatives, which
should be checked for mathematical completeness, but which are of minimal
practical importance. The result is a necessary condition for bifurcation..
These conditions are described in the following theorem. A sketch of the proof
can be found in [11], and a complete version is given in Appendix D.

######  Theorem 5.

Bifurcation occurs under the following condition

$\exists y_{n}\text{ s.t. }p(y_{n})>0\text{ and
}\det\left[I-\frac{1-\lambda}{\lambda}\frac{\partial^{2}\phi(y_{n})}{\partial
y_{n}^{2}}C_{X|y_{n}}\right]=0,$ (17)

where
$C_{X|y_{n}}:=\mathbb{E}\left[(X-y_{n})(X-y_{n})^{\mathrm{T}}|y_{n}\right]$.

###### Proof.

See Appendix D. ∎

In other words, there exist critical values for $\lambda$ that depend on the
data space itself and the choice of the Bregman divergence (through the
function $\phi$), such that bifurcation occurs when

$\frac{\lambda}{1-\lambda}=\frac{\partial^{2}\phi(y_{n})}{\partial
y_{n}^{2}}\bar{\nu}$ (18)

where $\bar{\nu}$ is the largest eigenvalue of $C_{X|y_{n}}$. That is to say
that an algorithmic implementation needs only as many codevectors as the
number of effective codevectors, which depends only on changes of the
temperature parameter below certain thresholds that depend on the dataset at
hand and the dissimilarity measure used. As shown in Alg. 1, we can detect the
bifurcation points by introducing perturbing pairs of codevectors at each
temperature level $\lambda$ (resp. $T$). In this way, the codevectors $\mu$
are doubled by inserting a perturbation of each $\mu_{i}$ in the set of
effective codevectors. The newly inserted codevectors will merge with their
pair if a critical temperature has not been reached and separate otherwise.

###  II-E Connection to Vector Quantization. Compression Rate and Error

It is apparent that problem (4) is an entropy-constrained generalization of a
soft-clustering method that, in the limit $\lambda\rightarrow 0$, converges to
a standard vector quantization (hard-clustering) problem. In fact, one can
easily verify that

$\displaystyle F_{\lambda}(\mu)+\lambda H(X)=(1-\lambda)D(\mu)-\lambda
H(\mu)+\lambda H(X)$ (19) $\displaystyle=\int
p(x)\sum_{i}p(\mu_{i}|x)\left[(1-\lambda)d(x,\mu_{i})-\lambda\log
p(\mu_{i}|x)\right]~{}dx,$

and, since the entropy term $H(X)$ does not depend on the optimization
parameters $\mu$, the clustering approach in (4) is equivalent to soft-
clustering with respect to a modified dissimilarity measure given by:

$d_{\lambda}(x,\mu_{i})=(1-\lambda)d(x,\mu_{i})-\lambda\log p(\mu_{i}|x)$ (20)

subject to the constraint $\sum_{i}p(\mu_{i}|x)=1$.

Therefore, the proposed method is a lossy compression method with
hierarchically decreasing loss as $\lambda\rightarrow 0$ and the number of
effective codevectors goes to infinity, i.e., $K\rightarrow\infty$. An
explicit expression of the error rate $F_{\lambda}(\mu^{*})$, for each
temperature level $\lambda$, as a function of $F_{0}(\mu^{*})=D(\mu^{*})$,
i.e., the error rate of a vector quantization algorithm, is hard to obtain as
it highly depends on the underlying distribution of the data space at hand
through the entropy term. However, an intuitive interpretation of the
hierarchy of solutions that is constructed by solving (4) for decreasing
values of $\lambda$ can be seen from the form of the conditional probabilities
in (6). That is, during the implementation of the algorithm, at every level
$\lambda$, a (soft-)Voronoi partition of the data space is computed with
respect to a scaled dissimilarity measure:

$\bar{d}_{\lambda}(x,\mu_{i})=\frac{1-\lambda}{\lambda}d(x,\mu_{i})=\frac{1}{T}d(x,\mu_{i})$
(21)

Thus, the algorithm perceives a scaled version of the data space at each level
$\lambda$, by focusing only to large dissimilarities within the data space
when the value of $\lambda$ is high, and progressively zooming in to perceive
more subtle dissimilarities as the value of $\lambda$ decreases. Therefore,
the error rate $F_{\lambda}(\mu^{*})$ can be roughly expressed as proportional
to $D(\mu^{*})$ and the term $\frac{1-\lambda}{\lambda}$ (inversely
proportional to the temperature level $T$), i.e.,
$F_{\lambda}(\mu^{*})\propto\frac{1-\lambda}{\lambda}D(\mu^{*})$.

Note that this is the worst case scenario, when the introduction of the
entropy term induces information loss across all regions of the data space. In
many cases, there are regions of the data space where higher compression rate
does not introduce information loss, or the information loss is significantly
lower than others. In this sense, one can view (4) as a risk-sensitive version
of soft-clustering, where an optimistic, or risk-seeking, approach is adopted.
Risk-seeking in this setting translates to searching for less complex
representations (with lower number of effective codevectors that induce higher
entropy) in the hope that more complex representations are not necessarily
needed. The parameter $\lambda$ then becomes a weight of risk-sensitivity.
More details regarding this interpretation can be found in [22]. In view of
the above, the error rate $F_{\lambda}(\mu^{*})$ can be roughly expressed as

$F_{\lambda}(\mu^{*})\leq\frac{1-\lambda}{\lambda}D(\mu^{*}),\
\lambda\in[0,1).$ (22)

###  II-F Algorithmic Implementation and Complexity

The progressive partitioning algorithm is shown in Algorithm 1. The
temperature parameter $\lambda_{t}$ is reduced using the geometric series
$\lambda_{t+1}=\gamma\lambda_{t}$, for $\gamma<1$. Regarding the stochastic
approximation stepsizes, simple time-based learning rates of the form
$\alpha_{n}=\nicefrac{{1}}{{a+bn}}$, $a,b>0$, have experimentally shown to be
sufficient for fast convergence. Convergence is checked with the condition
$\frac{1-\lambda}{\lambda}d_{\phi}(\mu_{i}^{n},\mu_{i}^{n-1})<\epsilon_{c}$
for a given threshold $\epsilon_{c}$. This condition becomes harder as the
value of $\lambda$ decreases. The stopping criteria $T_{stop}$ can include a
maximum number of codevectors $K_{max}$ allowed, a minimum temperature
$\lambda_{min}$ to be reached, a minimum distortion error $e_{target}$ to be
reached, a maximum number of iterations $i_{max}$, and so on.

Bifurcation, at $\lambda_{t}$, is detected by maintaining a pair
$\left\\{\mu_{j}+\delta,\mu_{j}-\delta\right\\}$ of perturbed codevectors for
each effective codevector $\mu_{j}$ generated by the algorithm at
$\lambda_{t-1}$, i.e. for $j=1\ldots,K_{i-1}$. Using arguments from
variational calculus (see Section II-D), it is easy to see that, upon
convegence, the perturbed codevectors will merge if a critical temperature has
not been reached, and will get separated otherwise. Therefore, the cardinality
of the model is at most doubled at every temperature level. These are the
effective codevectors discussed in Section II-D. Merging is detected by the
condition $\frac{1-\lambda}{\lambda}d_{\phi}(\mu_{j},\mu_{i})<\epsilon_{n}$,
where $\epsilon_{n}$ is a design parameter that acts as a regularization term
for the model that controls the number of effective codevectors. An additional
regularization mechanism is the detection of idle codevectors, which is
checked by the condition $\rho_{i}(n)<\epsilon_{r}$, where $\rho_{i}(n)$ can
be seen as an approximation of the probability $p(\mu_{i})$.

The complexity of Alg. 1 for a fixed temperature coefficient $\lambda_{t}$ is
$O(N_{c_{t}}(2K_{t})^{2}d)$, where $N_{c_{t}}$ is the number of stochastic
approximation iterations needed for convergence which corresponds to the
number of data samples observed, $K_{t}$ is the number of codevectors of the
model at temperature $\lambda_{t}$, and $d$ is the dimension of the input
vectors, i.e., $x\in\mathbb{R}^{d}$. Therefore, assuming a schedule
$\left\\{\lambda_{1}=\lambda_{max},\lambda_{2},\ldots,\lambda_{N_{\lambda}}=\lambda_{min}\right\\}$,
the time complexity for the training of Algorithm 1 becomes:

$\displaystyle O(N_{c}(2\bar{K})^{2}d)$

where $N_{c}=\max_{i}\left\\{N_{c_{t}}\right\\}$ is an upper bound on the
number of data samples observed until convergence at each temperature level,
and $\bar{K}=\sum_{i=1}^{N_{\lambda}}K_{t}$, with

$\displaystyle
N_{\lambda}\leq\bar{K}\leq\min\left\\{\sum_{n=0}^{N_{\lambda}-1}2^{n},\sum_{n=0}^{\log_{2}K_{max}}2^{n}\right\\}<N_{\lambda}K_{max}$

where the actual value of $\bar{K}$ depends on the bifurcations occurred as a
result of reaching critical temperatures and the effect of the regularization
mechanisms described above. Note that typically $N_{c}\ll N$ as a result of
the stochastic approximation algorithm, and $\bar{K}\ll N_{\lambda}K_{max}$ as
a result of the progressive nature of the algorithm. Prediction scales
linearly with $O(K_{N_{\lambda}}d)$, with $K_{N_{\lambda}}\leq K_{max}$.

Algorithm 1 Progressive Partitioning.

Select a Bregman divergence $d_{\phi}$

Set stopping criteria $T_{stop}$ (e.g., $K_{max}$, $\lambda_{min}$)

Set convergence parameters: $\gamma$, $\epsilon_{c}$, $\epsilon_{n}$,
$\epsilon_{r}$, $\delta$

Set stepsizes: $\left\\{\alpha_{n}\right\\}$

Initialize: $K=1$, $\lambda=1$,asdfasdfa $\left\\{\mu_{0}\right\\}$,
$p(\mu_{0})=1$, $\sigma(\mu_{0})=\mu_{0}p(\mu_{0})$

repeat

Perturb codebook:
$\left\\{\mu_{i}\right\\}\leftarrow\left\\{\mu_{i}+\delta\right\\}\bigcup\left\\{\mu_{i}-\delta\right\\}$

Update $K\leftarrow 2K$, $\left\\{p(\mu_{i})\right\\}$,
$\left\\{\sigma(\mu_{i})\leftarrow\mu_{i}p(\mu_{i})\right\\}$

$n\leftarrow 0$

repeat

Observe data point $x$

for $i=1,\ldots,K$ do

Update:  $\displaystyle p(\mu_{i}|x)$
$\displaystyle\leftarrow\frac{p(\mu_{i})e^{-\frac{1-\lambda}{\lambda}d_{\phi}(x,\mu_{i})}}{\sum_{i}p(\mu_{i})e^{-\frac{1-\lambda}{\lambda}d_{\phi}(x,\mu_{i})}}$
$\displaystyle p(\mu_{i})$ $\displaystyle\leftarrow
p(\mu_{i})+\alpha_{n}\left[p(\mu_{i}|x)-p(\mu_{i})\right]$
$\displaystyle\sigma(\mu_{i})$
$\displaystyle\leftarrow\sigma(\mu_{i})+\alpha_{n}\left[xp(\mu_{i}|x)-\sigma(\mu_{i})\right]$
$\displaystyle\mu_{i}$
$\displaystyle\leftarrow\frac{\sigma(\mu_{i})}{p(\mu_{i})}$

$n\leftarrow n+1$

end for

until Convergence:
$\frac{1-\lambda}{\lambda}d_{\phi}(\mu_{i}^{n},\mu_{i}^{n-1})<\epsilon_{c}$,
$\forall i$

Keep effective codevectors: asdf discard $\mu_{i}$ if
$\frac{1-\lambda}{\lambda}d_{\phi}(\mu_{j},\mu_{i})<\epsilon_{n}$, $\forall
i,j,i\neq j$

Remove idle codevectors: asdf discard $\mu_{i}$ if $p(\mu_{i})<\epsilon_{r}$,
$\forall i$

Update $K$, $\left\\{p(\mu_{i})\right\\}$, $\left\\{\sigma(\mu_{i})\right\\}$

Lower temperature: $\lambda\leftarrow\gamma\lambda$

until $T_{stop}$

##  III Learning with Local Models: Combined Partitioning and Function
Approximation

In this section, we investigate the problem of combined partitioning and
function approximation, which results in a learning approach where multiple
local models are trained, taking advantage of the differences in the
underlying probability distribution of the data space. As a consequence, this
approach can circumvent the use of overly complex learning models, reduce
time, memory, and computational complexity, and give insights to certain
properties of the data space [10].

In the general case, a function $f:S\rightarrow\mathcal{F}$ is to be
approximated given a set of observations $\left\\{(x_{n},f(x_{n}))\right\\}$
where $\left\\{x_{n}\right\\}$ are independent realizations of a random
variable $X\in S$, similar to Section II. One then seeks to find a partition
$\left\\{S_{i}\right\\}$ and a set of parameters
$\left\\{\theta_{i}\right\\}\in\Theta$ for some predefined learning models
$\left\\{\hat{f}_{i}(x,\theta_{i})\in\mathcal{F}\right\\}$ such that:

$\min_{\left\\{S_{i},\theta_{i}\right\\}}\
\mathbb{E}\left[\sum_{i}\mathds{1}_{\left[X\in
S_{i}\right]}d\left(f(X),\hat{f}_{i}(X,\theta_{i})\right)\right]$ (23)

where $d:\mathcal{F}\times\mathcal{F}\rightarrow[0,\infty)$ is a well-defined
convex metric with respect to the second argument.

To find a tractable solution to this problem, we decompose the two tasks of
progressive partitioning and function approximation. As described in Section
II, a partition $\left\\{S_{i}\right\\}_{i=1}^{K(\lambda)}$ of the space $S$
can be approximated online using a stochastic approximation algorithm that
solves (4) and yields the locations of a finite number of
$\mu_{\lambda}:=\left\\{\mu_{i}\right\\}_{i=1}^{K(\lambda)}$ codevectors, that
define the regions $S_{i}=\left\\{x\in
S:i=\operatorname*{arg\,min}_{j}~{}d_{\phi}(x,\mu_{j})\right\\}$,
$i=1,\ldots,K(\lambda)$. Given the partition $\left\\{S_{i}\right\\}$, we are
now in place to solve the following problem:

$\min_{\left\\{\theta_{i}\right\\}}\
\mathbb{E}\left[\sum_{i}\mathds{1}_{\left[X\in
S_{i}\right]}d\left(f(X),\hat{f}_{i}(X,\theta_{i})\right)\right]$ (24)

######  Remark 2.

Solving (24) decouples the two tasks of progressive partitioning and local
function approximation and yields a sub-optimal solution to the original
combined problem in (23). That being said, the use of Alg. 1 is a heuristic
method that offers (i) the crucial properties of progressive partitioning, and
(ii) a compressed representation of the data space $S$ such that each $S_{i}$
represents a region of $S$ where its underlying probability distribution
presents low variability (see Section II).

In the remaining section, we will study learning approaches to computationally
solve (24) in the general case of a differentiable (with respect to
$\theta_{i}$) learning model $\hat{f}_{i}(x,\theta_{i})$, and in the specific
case of using locally constant models, i.e., when
$\hat{f}(x,\theta_{i})=\theta_{i}\in\mathcal{F}$.

###  III-A Learning with Local Models

In this section, we assume a model $\hat{f}_{i}(x,\theta_{i})$
$\in\mathcal{F}$ that is differentiable with respect to a parameter vector
$\theta_{i}\in\Theta$, where $\Theta$ is a finite-dimensional vector space.
Given a finite partition set of parameters
$\left\\{S_{i}\right\\}_{i=1}^{K(\lambda)}$, for $K(\lambda)<\infty$, (24) is
decomposed to

$\min_{\theta_{i}}\ \mathbb{E}\left[\mathds{1}_{\left[X\in
S_{i}\right]}d\left(f(X),\hat{f}_{i}(X,\theta_{i})\right)\right],\
i=1,\ldots,K(\lambda).$ (25)

where $d:\mathcal{F}\times\mathcal{F}\rightarrow[0,\infty)$ is assumed a
metric that is differentiable and convex with respect to the second argument.

This is a stochastic optimization problem that can be solved using stochastic
approximation updates. In particular, since we have assumed that
$\hat{f}_{i}(x,\theta_{i})$ is differentiable with respect to $\theta_{i}$, we
can use stochastic gradient descent:

$\displaystyle\theta_{i}(n+1)=\theta_{i}(n)-\beta(n)\nabla_{\theta}d(f(x_{n}),\hat{f}_{i}(x_{n},\theta_{i}(n)))$
(26)
$\displaystyle\quad=\theta_{i}(n)-\beta(n)\\{\nabla_{\theta}\mathbb{E}\left[\mathds{1}_{\left[X\in
S_{i}\right]}d(f(x_{n}),\hat{f}_{i}(x_{n},\theta_{i}(n)))\right]$
$\displaystyle\quad\quad+\big{(}\nabla_{\theta}d(f(x_{n}),\hat{f}_{i}(x_{n},\theta_{i}(n)))$
$\displaystyle\quad\quad-\nabla_{\theta}\mathbb{E}\left[\mathds{1}_{\left[X\in
S_{i}\right]}d(f(x_{n}),\hat{f}_{i}(x_{n},\theta_{i}(n)))\right]\big{)}\\},\
x_{n}\in S_{i}$

Since we can control the observations for each model $f_{i}$ to belong to
$S_{i}$, it is easy to see that
$M_{n+1}:=\nabla_{\theta}d(f(x_{n}),\hat{f}_{i}(x_{n},\theta_{i}(n)))-\nabla_{\theta}\mathbb{E}\left[\mathds{1}_{\left[X\in
S_{i}\right]}d(f(x_{n}),\hat{f}_{i}(x_{n},\theta_{i}(n)))\right]$ is a
martingale difference sequence for an unbiased estimator
$\nabla_{\theta}d(f(x_{n}),\hat{f}_{i}(x,\theta_{i}))$, i.e., when the
condition $\mathbb{E}\left[\mathds{1}_{\left[X\in
S_{i}\right]}\nabla_{\theta}d(f(x_{n}),\hat{f}_{i}(x_{n},\theta_{i}(n)))\right]=\nabla_{\theta}\mathbb{E}\left[\mathds{1}_{\left[X\in
S_{i}\right]}d(f(x_{n}),\hat{f}_{i}(x_{n},\theta_{i}(n)))\right]$ holds.
Therefore, as an immediate result of Theorem 10 in Appendix B, the stochastic
approximation process (26) converges almost surely to a possibly path-
dependent invariant set of
$\dot{\theta}_{i}=\nabla_{\theta}\mathbb{E}\left[\mathds{1}_{\left[X\in
S_{i}\right]}\hat{f}_{i}(x,\theta_{i})\right]$, i.e., an asymptotically stable
local minimum of the objective function
$\mathbb{E}\left[\mathds{1}_{\left[X\in
S_{i}\right]}\hat{d}(f(x_{n}),f_{i}(x,\theta_{i}))\right]$.

So far, we have assumed that $\left\\{S_{i}\right\\}$ is fixed. However, we
are interested in a learning approach that approximates
$\left\\{S_{i}\right\\}$ and $\left\\{\hat{f}_{i}(x,\theta_{i})\right\\}$ at
the same time, and given the same observations
$\left\\{(x_{n},f(x_{n}))\right\\}$ which may be available one at a time (i.e,
no dataset is stored in memory a priori). This is possible because both
learning algorithms for $\left\\{S_{i}\right\\}$ and
$\left\\{\hat{f}_{i}(x,\theta_{i})\right\\}$ independently are stochastic
approximation algorithms. According to the theory of two-timescale stochastic
approximation, we can run both learning algorithms at the same time, but using
different stepsize profiles $\left\\{\alpha(n)\right\\}$ and
$\left\\{\beta(n)\right\\}$, such that
$\nicefrac{{\alpha(n)}}{{\beta(n)}}\rightarrow 0$. Intuitively, we create a
system of two dynamical system running in different “speed”, meaning that
second system, the one with stepsizes $\left\\{\beta(n)\right\\}$, is updated
fast enough that the first system, the one with stepsizes
$\left\\{\alpha(n)\right\\}$, can be seen as quasi-static with respect to the
second. The following theorem summarizes this result.

######  Theorem 6.

Let $\left\\{x_{n}\right\\}$ be a sequence of independent realizations of $X$,
and assume that $\mu_{i}(n)$ is a sequence updated using the stochastic
approximation algorithm in (12) with stepsizes $\left\\{\alpha(n)\right\\}$
satisfying $\sum_{n}\alpha(n)=\infty$, and $\sum_{n}\alpha^{2}(n)<\infty$.
Then, as long as $\left\\{\beta(n)\right\\}$ are designed such that
$\sum_{n}\beta(n)=\infty$, $\sum_{n}\beta^{2}(n)<\infty$, and
$\nicefrac{{\alpha(n)}}{{\beta(n)}}\rightarrow 0$, the asynchronous updates

$\theta_{i}(n+1)=\theta_{i}(n)-\beta(n)\nabla_{\theta}d(f(x_{n}),\hat{f}_{i}(x_{n},\theta_{i}(n))),$
(27)

for $i=\operatorname*{arg\,min}_{j}~{}d_{\phi}(x_{n},\mu_{j}(n))$ converges
almost surely to a locally asymptotically stable solution
$\left\\{\theta_{i}\right\\}$ of (25), as $n\rightarrow\infty$, for
$S_{i}=\\{x\in
S:i=\operatorname*{arg\,min}_{j}~{}d_{\phi}(x,\mu_{j}(\infty))\\}$, where
$\mu_{i}(\infty))$ is the asymptotically stable equilibrium of (12).

###### Proof.

See Appendix E. ∎

The algorithmic implementation is shown in Alg. 2 as an extension of Alg. 1.

Algorithm 2 Progressive Learning with Differentiable Models.

—–$//$—–

Set stepsizes: $\left\\{\alpha_{n}\right\\}$, $\left\\{\beta_{n}\right\\}$
s.t. $\nicefrac{{\alpha_{n}}}{{\beta_{n}}}\rightarrow 0$

Initialize: $\left\\{\mu_{0}\right\\}$, $\left\\{\theta_{0}\right\\}$

repeat

—–$//$—–

repeat

Observe data point $x$ & output $y$

for $i=1,\ldots,K$ do

Update:

—–$//$—–  $\displaystyle p(\mu_{i})$ $\displaystyle\leftarrow
p(\mu_{i})+\alpha_{n}\left[p(\mu_{i}|x)-p(\mu_{i})\right]$
$\displaystyle\sigma(\mu_{i})$
$\displaystyle\leftarrow\sigma(\mu_{i})+\alpha_{n}\left[xp(\mu_{i}|x)-\sigma(\mu_{i})\right]$
asfda
$\theta_{i}\leftarrow\theta_{i}-\beta_{n}\nabla_{\theta}d(f(x),\hat{f}_{i}(x,\theta_{i}))$

—–$//$—–

end for

until Convergence

—–$//$—–

until $T_{stop}$

###  III-B Case of Constant Local Models

In the special case when locally constant models are used, i.e., when
$\hat{f}(x,\theta_{i})=\theta_{i}\in\mathcal{F}$, two-timescale updates are
not required, and a simpler solution can be tracked. In particular, we can
augment the system (12) with

$\begin{cases}\sigma_{\theta_{i}}(n+1)&=\sigma_{\theta_{i}}(n)+\alpha(n)\left[x_{n}\hat{p}(\mu_{i}|x_{n})-\sigma_{\theta_{i}}(n)\right]\\\
\theta_{i}(n)&=\frac{\sigma_{\theta_{i}}(n)}{\rho_{i}(n)}\end{cases}$ (28)

Following the same arguments as in the proof of Theorem 3, it is easy to see
that $\theta_{i}(n)$ converge almost surely to
$\mathbb{E}\left[\mathds{1}_{\left[X\in S_{i}\right]}f(X)\right]$ as
$n\rightarrow\infty$ and $\lambda\rightarrow 0$. To see this, notice that as
$\lambda\rightarrow 0$, $p^{*}(x,\mu_{i})\rightarrow\mathds{1}_{\left[X\in
S_{i}\right]}$ and $p^{*}(\mu_{i})\rightarrow 1$. As a final note, this
approach is equivalent to a piece-wise constant approximation of $f(X)$. In
other words, this is a binning process where the size and location of the bins
depends on the underlying probability distribution of $X$, and the number of
bins progressively increases, resulting in a hierarchical approximation of
$f(X)$. The algorithmic implementation is shown in Alg. 3 as an extension of
Alg. 1.

Algorithm 3 Progressive Learning with Constant Models.

—–$//$—–

Initialize: $\left\\{\mu_{0}\right\\}$, $\left\\{f_{\mu_{0}}\right\\}$,
$\left\\{\sigma_{f}(\mu_{0})\right\\}$

repeat

—–$//$—–

repeat

Observe data point $x$ & output $y$

for $i=1,\ldots,K$ do

Update:

—–$//$—– asdfa
$\sigma_{f}(\mu_{i})\leftarrow\sigma(\mu_{i})+\alpha_{n}\left[yp(\mu_{i}|x)-\sigma(\mu_{i})\right]$

asdfa $\mu_{i}\leftarrow\frac{\sigma(\mu_{i})}{p(\mu_{i})}$,
$f_{\mu_{i}}\leftarrow\frac{\sigma_{f}(\mu_{i})}{p(\mu_{i})}$

end for

until Convergence

—–$//$—–

until $T_{stop}$

##  IV The Problem of Classification

In this section we focus on the binary classification problem. The results can
be extended to the general case (see, e.g., [23]). For the classification
problem, a pair of random variables $\left\\{X,c(X)\right\\}\in
S\times\left\\{0,1\right\\}$ defined in a probability space
$\left(\Omega,\mathcal{F},\mathbb{P}\right)$, is observed with $c(X)$
representing the class of $X$ and $S\subseteq\mathbb{R}^{d}$. The codebook is
represented by $\mu:=\left\\{\mu_{i}\right\\}_{i=1}^{K}$, $\mu_{i}\in ri(S)$,
and $c_{\mu}:=\left\\{c_{\mu_{i}}\right\\}_{i=1}^{K}$, such that
$c_{\mu_{i}}\in\left\\{0,1\right\\}$ represents the class of $\mu_{i}$ for all
$i\in\left\\{1,\ldots,K\right\\}$. A partition-based classifier is called
Bayes-optimal if it minimizes the classification error:

$\displaystyle\min_{\mu,c_{\mu}}~{}J_{B}(\mu,c_{\mu})$
$\displaystyle:=\pi_{1}\sum_{i:c_{\mu_{i}}=0}\mathbb{P}\left\\{X\in
S_{i}|c(X)=1\right\\}$ (29)
$\displaystyle\quad+\pi_{0}\sum_{i:c_{\mu_{i}}=1}\mathbb{P}\left\\{X\in
S_{i}|c(X)=0\right\\}$

where $S_{i}=\left\\{x\in
S:i=\operatorname*{arg\,min}\limits_{j}~{}d(x,\mu_{j})\right\\}$, and
$\pi_{i}:=\mathbb{P}\left[c=i\right]$.

In the remaining section, we study methods to solve the classification problem
based on the results of Sections II and III.

###  IV-A Classification as a Regression Problem with Constant Local Models

The classification problem (29) can be viewed as a special case of learning
with local models as in Section III. Here one seeks to find a partition
$\left\\{S_{i}\right\\}$ and a set of parameters
$\left\\{c_{i}\in\left\\{0,1\right\\}\right\\}$ such that:

$\min_{\left\\{S_{i},c_{i}\right\\}}\
\mathbb{E}\left[\sum_{i}\mathds{1}_{\left[X\in
S_{i}\right]}d\left(c(X),c_{i}\right)\right]$ (30)

where $d:=\mathds{1}_{\left[c\neq c_{\mu_{i}}\right]}$. Notice that since $d$
is not differentiable the results of Section III cannot be used directly.
However, numerous relaxation methods can be used to find a possibly sub-
optimal solution. A widely used approach is to relax the constraints on
$\left\\{c_{i}\in\left\\{0,1\right\\}\right\\}$ such that
$\left\\{c_{i}\in\left[0,1\right]\right\\}$. Then the updates (28) can be
directly used to estimate $\mathbb{E}\left[\mathds{1}_{\left[X\in
S_{i}\right]}c(X)\right]$. Then a projection mapping
$r:\left[0,1\right]\rightarrow\left\\{0,1\right\\}$, e.g.,
$r(c)=\mathds{1}_{\left[c<0.5\right]}$, can be used to return a solution to
the classification problem. Notice that this is equivalent to a majority-vote
rule inside each region $S_{i}$. This is a common approach that, at the limit
$\lambda\rightarrow 0$, when the updates (12), (13) result in a hard-
clustering approach with infinite number of clusters, yield a classification
rule that is strongly Bayes risk consistent, i.e., converges to the optimal
(Bayes) probability of error given in (29) (see, e.g., Ch. 21 in [23]).

###  IV-B Classification as Class-Conditioned Density Estimation

In a different approach, we can formulate the binary classification problem to
the minimization of $F$ in (4) with a modified average distortion measure
given by:

$\displaystyle D=\mathbb{E}\left[d^{b}(X,c,\mu,c_{\mu})\right]$

where $d^{b}(x,c,\mu_{i},c_{\mu_{i}})=\mathds{1}_{\left[x\in
S_{i}\right]}\mathds{1}_{\left[c\neq c_{\mu_{i}}\right]}$. However, because
$d^{b}$ is not differentiable, using similar principles as in the case of
Learning Vector Quantization (LVQ) [9], we can instead approximate the optimal
solution by using the distortion measure

$d^{l}(x,c_{x},\mu,c_{\mu})=\begin{cases}d(x,\mu),~{}c_{x}=c_{\mu}\\\
-d(x,\mu),~{}c_{x}\neq c_{\mu}\end{cases}$ (31)

Using similar arguments to Ch. 21 in [23], it can be shown that as
$\lambda\rightarrow 0$, the solution $(\mu,c_{\mu})$ to the above problem
equipped with a majority-vote classification rule is strongly Bayes risk
consistent.

However, we find useful to also explore a generative learning approach, using

$d^{c}(x,c_{x},\mu,c_{\mu})=\begin{cases}d(x,\mu),~{}c_{x}=c_{\mu}\\\
0,~{}c_{x}\neq c_{\mu}\end{cases}$ (32)

It is easy to see that this particular choice for the distortion measure
$d^{c}$ in (32) transforms the learning rule in (12) to

$\begin{cases}\rho_{i}(n+1)&=\rho_{i}(n)+\beta(n)\left[s_{i}\hat{p}(\mu_{i}|x_{n})-\rho_{i}(n)\right]\\\
\sigma_{i}(n+1)&=\sigma_{i}(n)+\beta(n)\left[s_{i}x_{n}\hat{p}(\mu_{i}|x_{n})-\sigma_{i}(n)\right]\end{cases}$
(33)

where $s_{i}:=\mathds{1}_{\left[c_{\mu_{i}}=c\right]}$. As a result, this is
equivalent to estimating strongly consistent class-conditional density
estimators:

$\hat{p}(x|c=j)\rightarrow\pi_{j}p(x|c=j),\ a.s.$ (34)

and the following theorem holds:

######  Theorem 7.

In the limit $\lambda\rightarrow 0$, and as the number of observed samples
$\left\\{x_{n}\right\\}$ goes to infinity, i.e., $n\rightarrow\infty$, the
learning algorithm based on (33), (13), results in strongly consistent class-
conditional density estimators $\hat{p}(x|c=j)$ that construct a Bayes risk
consistent classifier with the classification rule

$\hat{c}=\operatorname*{arg\,max}_{j}\hat{\pi}_{j}\hat{p}(x|c=j),\ j=1,2$ (35)

where $\hat{\pi}_{j}=\frac{\sum_{n}\mathds{1}_{\left[c_{n}=j\right]}}{n}$

###### Proof.

See Appendix F. ∎

As a final note, an easy-to-implement nearest-neighbor rule classification
rule:

$\hat{c}(x)=c_{\mu_{h^{*}}}$ (36)

where
$h^{*}=\operatorname*{arg\,max}\limits_{\tau=1,\ldots,K}~{}p(\mu_{\tau}|x),~{}h\in\left\\{1,\ldots,K\right\\}$,
yields a classification error $\hat{J}_{B}^{*}$ with tight upper bound with
respect to the Bayes-optimal $J_{B}^{*}$, i.e.,
$J_{B}^{*}\leq\hat{J}_{B}^{*}\leq 2J_{B}^{*}$ (see, e.g., [23]). The
algorithmic implementation is shown in Alg. 4 as an extension of Alg. 1.

Algorithm 4 Progressive Classification via Class-Conditional Density
Estimation.

—–$//$—–

Initialize: $\left\\{\mu_{0}\right\\}$ & $\left\\{c_{\mu_{0}}\right\\}$,

repeat

—–$//$—–

repeat

Observe data point $x$ & class label $c$

if $\nexists\mu_{i}$ s.t. $c_{\mu_{i}}=c$ then

Insert:
$\left\\{\mu_{i}\right\\}\leftarrow\left\\{\mu_{i}\right\\}\bigcup\left\\{x\right\\}$

Insert:
$\left\\{c_{\mu_{i}}\right\\}\leftarrow\left\\{c_{\mu_{i}}\right\\}\bigcup\left\\{c\right\\}$

end if

for $i=1,\ldots,K$ do

Compute membership $s_{i}=\mathds{1}_{\left[c_{\mu_{i}}=c\right]}$

Update:

—–$//$—– asdfd $p(\mu_{i})\leftarrow
p(\mu_{i})+\alpha_{n}\big{[}$$s_{i}$$p(\mu_{i}|x)-p(\mu_{i})\big{]}$ asdfd
$\sigma(\mu_{i})\leftarrow\sigma(\mu_{i})+\alpha_{n}\big{[}$$s_{i}$$xp(\mu_{i}|x)-\sigma(\mu_{i})\big{]}$

—–$//$—–

end for

until Convergence

—–$//$—–

until $T_{stop}$

##  V Hierarchical Learning in Multiple Resolutions

In this section, we extend the progressive partitioning algorithm (Alg. 1) of
Section II, by imposing a tree structure in the construction of the regions
$\left\\{S_{i}\right\\}$. The results of Sections III, IV, are extended
naturally through their immediate dependence on the partition
$\left\\{S_{i}\right\\}$. The key idea of the progressive construction of the
tree structure is as follows. Given a value for the temperature coefficient
$\lambda_{t}$, Algorithm 1, as presented in Section II, yields a sequence of
partitions $\left\\{S_{i}\right\\}_{i=1}^{K(\lambda^{t})}$. If at
$\lambda^{t}$, a user-defined splitting criterion is met, the partition
$\left\\{S_{i}\right\\}_{i=1}^{K(\lambda^{t})}$ is fixed, and Algorithm 1 is
applied independently to each region $S_{i}$ to create
$\left\\{\left\\{S_{ij}\right\\}_{j=1}^{K_{i}(\lambda^{t_{i}})}\right\\}_{i=1}^{K(\lambda^{t})}$,
such that $\left\\{S_{ij}\right\\}_{j=1}^{K_{i}(\lambda^{t_{i}})}$ form a
partition of $S_{i}$. This is depicted in Fig. 1(c). For each parent set
$S_{i}$, the number of children sub-sets $K_{i}(\lambda^{t_{i}})$, may be
different, depending on the properties of $S_{i}$. The same holds for the
stopping values $\lambda^{t_{i}}$. The splitting criterion can involve terms
such as a minimum value of $\lambda_{min}$ reached, a maximum number of
$K_{max}$ codevectors reached, or a minimum percentage of improvement in
accuracy or distortion reduction for every temperature step is reached.

This structural constraint reduces the time complexity of the algorithm from
$O(K^{2})$ to $O(k^{2}+\log_{k}K)$, where $K$ here represents the total number
of sets $\left\\{S_{i}\right\\}_{i=1}^{K}$, and $k$ represents the number of
children sub-sets for each parent set (assumed equal for every parent set)
[24]. In addition, as we will show, this tree structure offers an inherent
regularization mechanism in classification applications (Section IV). Finally,
since the resulting structure is a non-binary tree-structure, we are able to
control the number of layers of the tree-structured partition of the data
space, without sacrificing the performance of the learning algorithm, i.e., a
finite tree depth is sufficient for convergence [25, 26]. Therefore, we can
match the number of tree layers to the number of resolutions in a multi-
resolution data representation. This will allow for training each layer of the
tree with progressively finer resolution of the data representation, which
defines a hierarchical and progressive learning approach that further reduces
the complexity of the algorithm, while inheriting potential benefits from the
feature extraction process of the multi-resolution analysis. As we will show,
in the case when group-convolutional wavelet transform is used to create the
multi-resolution data representation, this architecture shares similar
properties to a deep neural network architecture [27].

###  V-A Tree-Structured Progressive Partitioning

A tree-structured partition $\Sigma_{\Delta}:=\left\\{S_{\nu_{i}}\right\\}$ is
defined by a set of regions $S_{\nu_{i}}\in S$, each represented by a tree
node $\nu_{i}$, arranged in a tree structure $\Delta$ with a single root node
$\nu_{0}$ such that $S_{\nu_{0}}=S$. The tree structure $\Delta$ is a special
case of a connected, acyclic directed graph, where each node has a single
parent node (except for the root node) and an arbitrary number of children
nodes, that is, $\Delta$ is not restricted to be a binary tree. The set
$C(\nu_{i})$ represents the nodes $\left\\{\nu_{j}\right\\}$ that are children
of $\nu_{i}$, while the set $P(\nu_{j})$ represents the node $\nu_{i}$ for
which $\nu_{j}\in C(\nu_{i})$. The level $l\geq 0$ of a node
$\nu_{h}\in\Delta$ is the length of the path
$\left\\{\nu_{0},\ldots,\nu_{i},\nu_{j},\ldots,\nu_{h}\right\\}$ leading from
the root node $\nu_{0}$ to $\nu_{h}$ such that $\nu_{j}\in C(\nu_{i})$. The
terminal nodes $\tilde{\nu}:=\left\\{\nu_{i}:C(\nu_{i})=\emptyset\right\\}$
are called leaves, and the union of their associated sets will be denoted
$\tilde{S}:=\left\\{\tilde{S}_{j}\right\\}$, where $|\tilde{S}|=\tilde{K}$ is
the number of leaf sets that create a partition of $S$, and
$\tilde{l}:=\max\left\\{l:\nu_{i}^{(l)}\in\tilde{\nu}\right\\}<\infty$ will
denote the maximum depth of the tree.

$\Sigma_{\Delta}$ defines a hierarchical partitioning scheme for the domain
$S$, such that for every node $\nu_{i}\in\Delta$ associated with the region
$S_{\nu_{i}}$, its children nodes $\left\\{\nu_{j}\in C(\nu_{i})\right\\}$ are
associated with the regions $\left\\{S_{\nu_{j}}\right\\}$ that form a
partition of $S_{\nu_{i}}$. We will use the unique paths from the root node as
identification label for each node, i.e., $\nu_{j}=0\ldots ij$ such that
$\nu_{0}=0$, $C(0)=\left\\{0i\right\\}$, $C(0i)=\left\\{0ij\right\\}$, and so
on. As such, Algorithm 1 can be used recursively to construct a tree-
structured partition $\Sigma_{\Delta}$ as follows: Start with node $\nu_{0}=0$
as the only leaf node. Using observations $\left\\{x_{n}\right\\}$
(realizations of $X\in S$), apply Algorithm 1 until a partition
$\left\\{S_{0j}\right\\}$ of $S_{0}=S$ is constructed. Then starting with
$w=0$ and for every observation $x_{n}$, iterate the process

$\text{repeat }w\leftarrow w^{\prime}\in C(w)\text{ such that }x_{n}\in
S_{w^{\prime}}\text{, until }C(w)=\emptyset,$ (37)

and apply one stochastic approximation update of Algorithm 1 in $S_{w}$. This
asynchronous process can continue until the convergence of all applications of
Alg. 1, when a finite-depth tree-structured partition $\Sigma_{\Delta}$ is
constructed such that for every node $w\in\Delta$ with children nodes
$\left\\{wj\right\\}\in C(w)$, the regions $\left\\{S_{wj}\right\\}$ form a
partition of $S_{w}$. This process is illustrated in Algorithm 5 and its
asymptotic behavior is given by the following theorem:

######  Theorem 8.

Let $\Sigma_{\Delta}$ be a finite-depth tree-structured partitioning scheme
created by Alg. 5 using realizations $\left\\{x_{n}\right\\}$ of a random
variable $X\in S$. If the leaf nodes are updated at the limit
$\lambda\rightarrow 0$, and $n\rightarrow\infty$, then $\Sigma_{\Delta}$
yields a consistent density estimator of $X$, with
$\hat{p}(x)=\frac{\sum_{i}\mathds{1}_{\left[x\in\tilde{S}_{i}\right]}}{nVol(\tilde{S}_{i})}$,
where $\tilde{S}_{i}$ is a leaf node given by the iterative process (37).

###### Proof.

It follows directly by the application of Theorem 4 to each region
$\tilde{S}_{j}$, where $\left\\{\tilde{S}_{j}\right\\}$ are the leaf nodes of
$\Sigma_{\Delta}$ that form a partition of $S$. ∎

######  Remark 3.

We have shown the tree-structured extension of Alg. 1, as well as its
asymptotic behavior for finite tree depth. It is straightforward to show that
similar results hold for the tree-structured extension of Algorithms 2, 3, and
4, regarding the regression and classification problems discussed in Sections
III, and IV.

Notice that a finite tree depth is sufficient for convergence to a consistent
density estimator. This follows from the fact that the children of each tree
node $\nu_{i}$, representing a region $S_{\nu_{i}}\in S$, are the output of
the progressive construction of a partition of $S_{\nu_{i}}$, based on Alg. 1.
This result will be used in Section V-B to build a multi-resolution extension
of Algorithm 5.

The time complexity of the tree-structured algorithm is significantly reduced.
Let $K_{max}$ be the total number of codevectors allowed. Then Alg. 1 has a
worst-case complexity $O(N_{c}(2\bar{K})^{2}d)$, for training, where
$\bar{K}=\sum_{n=0}^{\log_{2}K_{max}}2^{n}$, while testing requires
$O(K_{max}d)$ (see Section II-F for details on the parameters). In a tree-
structured partition $\Sigma_{\Delta}$ of depth $\tilde{l}$, assuming that the
number of children $k=|C(\nu_{i})|$ of each node $\nu_{i}$ is $k$ is the same,
and that each region is represented by roughly the same number of observations
$N_{c}^{\tilde{l}}$, we get $k=(K_{max})^{\nicefrac{{1}}{{\tilde{l}}}}$, and
$N_{c}^{\tilde{l}}=\nicefrac{{N_{c}}}{{k}}$. Then training requires in the
worst-case:

$\displaystyle
O\left(\frac{k^{\tilde{l}}-1}{k(k-1)}N_{c}(2\bar{k})^{2}d)\right)$

where
$\bar{k}=\sum_{n=0}^{\log_{2}k}2^{n}=\sum_{n=0}^{\nicefrac{{1}}{{\tilde{l}}}\log_{2}K_{max}}2^{n}$.
Prediction requires a forward pass of the tree, i.e., it scales with
$O(k\log_{k}K_{max}d)$.

In addition, we note that Alg. 5 updates the partition
$\left\\{S_{wj}\right\\}$ of each node $w$ asynchronously. As a result,
depending on the underlying probability density of the random variable $X\in
S$, some nodes will be visited more often than others, which will result in
some branches of the tree growing faster than others, inducing a variable-rate
code that frequently outperforms fixed-rate, full-search techniques with the
same average number of bits per sample. Alternatively, when learning offline
using a dataset, all nodes can be trained using parallel processes, which can
be utilized by multi-core computational units.

In the classification problem, an additional regularization mechanism can be
added in the approach described in Section IV-B (Alg. 4). Specifically, when a
partition $\left\\{S_{wj}\right\\}_{j=1}^{K_{w}}$ of a node $w$ is fixed, the
node $w$ can check the condition $c_{\mu_{wi}}=c_{\mu_{wj}},\forall i,j$,
which means that the partition $\left\\{S_{wj}\right\\}_{j=1}^{K_{w}}$ is
using $K_{w}$ codevectors, all of which correspond to the same class. In this
case, node $w$ is assigned a single codevector, and is not further split by
the algorithm. This phenomenon is illustrated in Fig. 7 in Section VI.

Finally, we note that the termination criteria of the iterations of Alg. 1 in
Alg. 5 in each layer of the tree are important design parameters. These can
include a maximum number of codevectors for each partition, a minimum
temperature $\lambda_{min}^{l}$ in each tree layer $l$, a maximum number of
iterations, and so on. These termination criteria characterize the splitting
criteria as well, i.e., when to stop growing the set of effective codevectors
in a partition of layer $l$, and continue to split the partition in the next
layer $l+1$.

Algorithm 5 Tree-Structured Progressive Partitioning

Initialize root node $\nu_{0}$ s.t. $S_{\nu_{0}}=S$

repeat

Observe data point $x\in S$

Find leaf node to update:

Set $w=\nu_{0}$

while $C(w)\neq\emptyset$ do

$w\leftarrow v\in C(w)$ such that $x\in S_{v}$

end while

Update partition $\left\\{S_{wj}\right\\}$ of $S_{w}$ using $x$ and Alg. 1

if Alg. 1 in $S_{w}$ terminates then

Split node $w$: $C(w)\leftarrow\left\\{S_{wj}\right\\}$

end if

until Stopping criterion

###  V-B Multi-Resolution Extension

So far we have modeled the observations as realizations of a random variable
$X\in S\subseteq\mathbb{R}^{d}$. In general, $X$ can be itself a measurable
signal $X(t):\mathbb{R}^{n}\rightarrow S$ with finite energy, i.e., $X(t)\in
S\subseteq L^{2}(\mathbb{R}^{d})$. We will denote the original space $S$ as
$S^{0}$. A multi-resolution representation of the signal $X(t)$ consists of a
sequence of projections of $X(t)$ on subspaces $\left\\{S^{j}\right\\}$ such
that $S^{j}\subset S^{j-1}$, $\forall j\in\mathbb{N}$, and
$\cup_{j=0}^{\infty}S^{j}$ is dense in $S^{0}$ with
$\cap_{j=0}^{\infty}S^{j}=\left\\{0\right\\}$. There are numerous methods to
construct subspaces $\left\\{S^{j}\right\\}$ with these properties, from the
classical wavelet transform [28] to different dictionary learning approaches
[29, 30]. An approach using group-convolutional wavelet decomposition will be
presented in Section V-C.

We denote by $X^{r}\in S^{r}$ the projection of $X=X^{0}$ to the subspace
$S^{r}$. Given a multi-resolution representation of $X$ with subspaces
$\left\\{S^{0},S^{1},S^{2},\ldots,S^{\tilde{l}}\right\\}$, we can extend Alg.
5 presented in Section V-A such that $X^{\tilde{l}-r}\in S^{\tilde{l}-r}$ is
used to train the nodes of the tree at level $r$. This idea matches the
intuition of using higher-resolution representation of $X$ for deeper layers
of the tree. It was first introduced in [31] and [32], and constitutes a
hierarchical multi-resolution learning scheme. The algorithmic implementation
is straightforward given Alg. 5, and is given in Alg. 6.

Regarding the asymptotic behavior of the multi-resolution extension, the
results of Theorem 8 hold. To see that, notice that since $S^{l}\subset
S^{0}$, for $l>0$, it follows that $X^{l}\in S^{0}$ as well, and, as a result,
Alg. 6 essentially creates a tree-structured partition $\Sigma_{\Delta}$ of
$S^{0}=S$, with the leaf nodes trained with $X^{0}\in S^{0}=S$, such that the
following holds.

######  Theorem 9.

Let $\Sigma_{\Delta}$ be a tree-structured partitioning scheme of depth
$\tilde{l}$ created by Alg. 6 using the multi-resolution representation
$\left(x_{n}^{0}=x_{n},x_{n}^{1},\ldots,x_{n}^{\tilde{l}}\right)\in\left(S^{0}=S,S^{1},\ldots,S^{\tilde{l}}\right)$
of realizations of a random variable $X\in S$. If the leaf nodes are updated
at the limit $\lambda\rightarrow 0$, and $n\rightarrow\infty$, then
$\Sigma_{\Delta}$ yields a consistent density estimator of $X\in S$, with
$\hat{p}(x)=\frac{\sum_{i}\mathds{1}_{\left[x\in\tilde{S}_{i}\right]}}{nVol(\tilde{S}_{i})}$,
where $\tilde{S}_{i}$ is a leaf node given by the iterative process (37).

###### Proof.

It follows directly by the application of Theorem 4 to each region
$\tilde{S}_{j}$, where $\left\\{\tilde{S}_{j}\right\\}$ are the leaf nodes of
$\Sigma_{\Delta}$ that form a partition of $S$. ∎

Algorithm 6 Multi-Resolution Progressive Partitioning

Initialize root node $\nu_{0}$ s.t. $S_{\nu_{0}}=S^{\tilde{l}}$

repeat

Observe data point $x^{\tilde{l}}\in S^{\tilde{l}}$

Find leaf node to update:

Set $w=\nu_{0}$

Set resolution $l=\tilde{l}$

while $C(w)\neq\emptyset$ do

$w\leftarrow v\in C(w)$ such that $x^{l}\in S_{v}^{l}$

$l\leftarrow l-1$

end while

Update partition $\left\\{S_{wj}^{l}\right\\}$ of $S_{w}^{l}$ using $x$ and
Alg. 1

if Alg. 1 in $S_{w}^{l}$ terminates and $l>0$ then

Split node $w$: $C(w)\leftarrow\left\\{S_{wj}^{l}\right\\}$

end if

until Stopping criterion

(a) DCNN.

(b) SCN.

(c) Proposed Architecture.

Figure 2: Block-diagram of the proposed hierarchical architecture using multi-
resolution features from the wavelet scattering transform compared to Deep
Convolutional Neural Networks (DCNN) and Scattering Convolutional Networks
(SCN). The feed-forward arrows represent a cascade of convolution, rectifying,
and downsampling operations.

###  V-C Building Group-Invariant Multi-Resolution Representations

There are numerous methods to construct subspaces $\left\\{S^{j}\right\\}$
with the properties mentioned in Section V-A. In this section we briefly
mention a particular approach based on group-convolutional wavelet
decomposition that aligns with the principles of the scattering transform,
first introduced in [19]. This is an unsupervised method that constructs a
hierarchy of features based on group-convolutions, that preserve local
invariance with respect to a certain class of Lie groups, such as translation,
rotation, and deformation, an important property in many learning applications
[27].

We start with the standard wavelet transform $\left\\{W_{l}X(t)\right\\}_{l}$
of a signal $X(t)\in S\subseteq L^{2}(\mathbb{R}^{d})$ as a basis [28]. Here,
$\left\\{W_{l}X(t)\right\\}_{l}\in S^{l}$ represents the signal $X(t)$ at
resolution $2^{l}$. The computation of the multi-resolution wavelet
representation of a signal consists of successive operations of a linear
convolution operator, followed by a downsampling step [27]. As a result, the
wavelet transform, is stable to small deformations [19]. In addition, the
wavelet transform is translation covariant (or equivariant), that is it
commutes with the Lie group of operators
$\left\\{T_{c}\right\\}_{c\in\mathbb{R}}$ such that $T_{c}X(t)=X(t-c)$, i.e.,
$W_{j}(T_{c}X)=T_{c}W_{j}(X)$. We note that, in the control theory and signal
processing communities, convolutions are associated with systems described by
the term ’linear time-invariant’. To avoid confusion, with the terminology
used here, these systems are considered linear covariant (or equivariant)
operators with respect to translation in time.

To induce local invariance (up to a scale $2^{J}$ for some $J>0$) with respect
to translation, it has been shown that it is sufficient to cascade the wavelet
transform with a non-linear operation $\rho W_{j}X=\|W_{j}X\|_{1}$, and a
locally averaging integral operation which can be modeled as a convolution
with a low-pass filter localized in a spatial window scaled at $2^{J}$ [19].
This is called a scattering transform and its implementation is based on a
complex-valued convolutional neural network whose filters are fixed wavelets
and $\rho$ is a complex modulus operator as described above [33]. This
structure is similar to deep convolutional neural networks [2], where
successive operations of a linear convolutional operator, a nonlinear mapping
(often a rectifying function, e.g., ReLu), and a down-sampling step (e.g.,
max-pooling), are used to produce the input for the next stage of the
architecture [19, 27, 34]. We illustrate this in Fig. 2, where Alg. 6 combined
with the hierarchical representation of a scattering transform is compared to
a Deep Convolutional Network [2] and a Scattering Convolutional Network [19].

As a final note, the translation invariance properties discussed above can be
generalized to the action of arbitrary compact Lie groups [35]. In particular,
let $G$ be a compact Lie group and $L^{2}(G)$ be the space of measurable
functions $f(r)$ such that $\|f\|^{2}=\int_{G}|f(r)|^{2}dr<\infty$, where $dr$
is the Haar measure of $G$. The left action of $g\in G$ on $f\in L^{2}(G)$ is
defined by $L_{g}f(r)=f(g^{-1}r)$. As a special case, the action of the
translation group $T_{c}f(t)=f(t-c)$ translates the function $f$ to the right
by $c$, with $g^{-1}=-c$ translating the argument of $f$ to the left by $c$.
Similar to the usual convolution $(f\ast
h)(x)=\int_{-\infty}^{\infty}f(u)h(x-u)du$ that defines a linear translation
covariant operator, convolutions on a group appear naturally as linear
operators covariant to the action of a group:

$(f\ast h)(x)=\int_{G}f(g)h(g^{-1}x)dr$ (38)

where $dr$ is the Haar measure of $G$. As a result, an invariant
representation relatively to the action of a compact Lie group, can be
computed by averaging over covariant representations created by group
convolution with appropriately defined wavelets, similar to the methodology
explained above.

##  VI Experimental Evaluation and Discussion

We illustrate the properties and evaluate the performance of the proposed
learning algorithm in clustering, classification, and regression problems.

In Fig. 3, the evolution of the progressive partitioning algorithm (Alg. 1)
studied in Section II is depicted, in an unsupervised learning (clustering)
problem. To better illustrate the properties of the approach, the data samples
were sampled from a mixture of 2D Gaussian distributions. The temperature
level (we use $T$ instead of $\lambda$ to stress the connection to the
temperature level in annealing optimization), the average distortion of the
model, the number of codevectors (neurons) used, the number of observations
(data samples) used for convergence, as well as the overall time, are shown.
This process showcases the performance-complexity trade-off described in
Section II. In Fig. 4 and 5, the tree-structured progressive partitioning
algorithm of Section V-A is compared against Alg. 1 in the same problem as in
Fig. 3. Notice that the time complexity of the algorithm is drastically
reduced. Additional properties regarding the construction of tree-structured
partitions are discussed in Section VI-B.

(a) Evolution of the algorithm in the data space.

(b) Performance curves.

Figure 3: Performance curves and data space evolution of Algorithm 1 applied
to a clustering problem with underlying Gaussian distributions.

(a) Evolution of the algorithm in the data space.

(b) Performance curves.

Figure 4: Performance curves and data space evolution of the tree-structured
approach (two layers) applied to a clustering problem with underlying Gaussian
distributions.

(a) Evolution of the algorithm in the data space.

(b) Performance curves.

Figure 5: Performance curves and data space evolution of the tree-structured
approach (three layers) applied to a clustering problem with underlying
Gaussian distributions.

Similarly, Fig. 6 shows the evolution of the learning model for a 2D
classification problem with class-conditional distributions given by a mixture
of 2D Gaussians. These results correspond to the approach explained in Section
IV-B, using Alg. 4. In addition to the apparent accuracy-complexity trade-off,
we make use of this classification problem to showcase the difference of using
the tree-structured approach of Section V-A, in terms of computational
complexity. In Fig. 7, we illustrate the evolution of the tree-structured
approach (Alg. 5). There are two notable comments on the behavior of this
approach compared to the original. First, the time complexity is considerably
improved (see Section V), and this results in a drastic difference in the
running time of the learning algorithm. Secondly, the number of codevectors
used is drastically reduced, as well. This is due to the regularization
mechanism described in Section V-A: when a partition
$\left\\{S_{wj}\right\\}_{j=1}^{K_{w}}$ of a node $w$ is fixed, the node $w$
can check the condition $c_{\mu_{wi}}=c_{\mu_{wj}},\forall i,j$, which means
that the partition $\left\\{S_{wj}\right\\}_{j=1}^{K_{w}}$ is using $K_{w}$
codevectors, all of which correspond to the same class. In this case, node $w$
is assigned only a single codevector, and is not further split by the
algorithm. Notice, that, in this way, the codevectors created by the tree-
structured algorithm tend to exist in the boundaries of the Bayes decision
surface, instead of populating areas where the decision surface does not
fluctuate at all.

(a) Evolution of the algorithm in the data space.

(b) Performance curves.

Figure 6: Performance curves and data space evolution of the proposed
algorithm applied to a classification problem with underlying Gaussian
distributions.

(a) Evolution of the algorithm in the data space.

(b) Performance curves.

Figure 7: Performance curves and data space evolution of the proposed tree-
structured algorithm (three layers) applied to a classification problem with
underlying Gaussian distributions.

The effect of using multiple resolutions as described in Section V-B, is
depicted in Fig. 8 for the same problem as in Fig. 6 and 7. For better
visualization, we assume that the low-resolution features, with respect to
which the first layer of the tree is computed, are the projections of the two-
dimensional data in an one-dimensional space (line). The second layer of the
tree is trained using the high-resolution features, i.e., the full knowledge
of both coordinates of the data. Notice that, as expected from Theorem 9, this
process will converge to a consistent learning algorithm, as long as the
multi-resolution representation used complies with the properties mentioned in
Section V-B, and the last layer of the tree uses the full knowledge of the
input data.

(a) Convergence of first layer with low-resolution features.

(b) Convergence of second layer with high-resolution features.

(c) Performance curve.

Figure 8: Performance curves and data space evolution of the proposed multi-
resolution algorithm (two layers) applied to a classification problem with
underlying Gaussian distributions.

Finally, in Fig. 9 and 10, we test the proposed methodology in two regression
problems, where one- and two-dimensional functions are hierarchically
approximated using the piece-wise constant approximation algorithm of Section
III-B, and the tree-structured approach of Section V.

(a) Evolution of the algorithm in the data space.

(b) Performance curves.

Figure 9: Performance curves and data space evolution of the proposed tree-
structured algorithm (four layers) applied to a piece-wise constant function
approximation problem in 1D.

(a) Evolution of the algorithm in the data space (original function on the
right).

(b) Performance curves.

Figure 10: Performance curves and data space evolution of the proposed tree-
structured algorithm (three layers) applied to a piece-wise constant function
approximation problem in 2D.

###  VI-A Source Code and Reproducibility

The open-source code is publicly available at
https://github.com/MavridisChristos/OnlineDeterministicAnnealing.

###  VI-B Tree-Structured Partition, Localization, and Explainability in
Machine Learning

Another advantage of using a tree-structured learning module is the
localization properties which allow for an understanding of the input space,
in accordance to the principles of the recently intoduced class of explainable
learning models [36]. The Voronoi regions shrink geometrically, and allow for
the use of local models, which is especially important in high-dimensional
spaces. Unlike most learning models, it is possible to locate the area of the
data space that presents the highest error rate and selectively split it by
using local ODA. This process can be iterated until the desired error rate (or
average distortion) is achieved. When using a training dataset for
classification, it is often possible to force accuracy of up to $100\%$ on the
training dataset. This is similar to an over-fitted classification and
regression tree (CART) [37]. However, over-fitting on the training dataset
often adversely affects the generalization properties of the model, the
performance on the testing dataset, and the robustness against adversarial
attacks. Therefore, the progressive process of ODA becomes important in
establishing a robust way to control the trade-off between performance and
complexity, before you reach that limit. Finally, an important question in
tree-structured learning models is the question of which cell to split next.
An exhaustive search in the entire tree to find the node that presents the
largest error rate is possible but is often not desired due to the large
computational overhead. This is automatically answered by the multi-resolution
ODA algorithm (Alg. 6) as it asynchronously updates all cells depending on the
sequence of the online observations. As a result, the regions of the data
space that are more densely populated with data samples are trained first,
which results in a higher percentage of performance increase per cell split.
We stress that this property makes the proposed algorithm completely dataset-
agnostic, in the sense that it does not require the knowledge of a training
dataset a priori, but instead operates completely online, i.e., using one
observation at a time to update its knowledge base.

##  VII Conclusion

We introduced a hierarchical learning algorithm to gradually approximate a
solution to a data-driven optimization problem in the context of autonomous
decision-making systems, especially under limitations on time and
computational resources. The learning architecture simulates an annealing
process and defines a heuristic method to progressively construct a tree-
structured partition of a possibly multi-resolution data space, which can be
used in conjunction with general learning algorithms to train local models.
The structured partitioning of the input space provides explainability, and
makes the learning architecture a suitable candidate for transfer learning
applications. Finally, the online gradient-free training rule based on
stochastic approximation, can be viewed a discrete-time dynamical learning
system, and used for inference, control, and reinforcement learning
applications.

## References

  * [1] K. P. Bennett and E. Parrado-Hernández, “The interplay of optimization and machine learning research,” _The Journal of Machine Learning Research_ , vol. 7, pp. 1265–1281, 2006.
  * [2] Y. LeCun, Y. Bengio, and G. Hinton, “Deep learning,” _nature_ , vol. 521, no. 7553, pp. 436–444, 2015.
  * [3] G. E. Hinton, S. Osindero, and Y.-W. Teh, “A fast learning algorithm for deep belief nets,” _Neural computation_ , vol. 18, no. 7, pp. 1527–1554, 2006\.
  * [4] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” in _Advances in neural information processing systems_ , 2012, pp. 1097–1105.
  * [5] H. Lee, R. Grosse, R. Ranganath, and A. Y. Ng, “Convolutional deep belief networks for scalable unsupervised learning of hierarchical representations,” in _Proceedings of the 26th annual international conference on machine learning_ , 2009, pp. 609–616.
  * [6] N. C. Thompson, K. Greenewald, K. Lee, and G. F. Manso, “The computational limits of deep learning,” _arXiv preprint arXiv:2007.05558_ , 2020.
  * [7] E. Strubell, A. Ganesh, and A. McCallum, “Energy and policy considerations for deep learning in nlp,” _arXiv preprint arXiv:1906.02243_ , 2019.
  * [8] M. Biehl, B. Hammer, and T. Villmann, “Prototype-based models in machine learning,” _Wiley Interdisciplinary Reviews: Cognitive Science_ , vol. 7, no. 2, pp. 92–111, 2016.
  * [9] T. Kohonen, _Learning Vector Quantization_. Berlin, Heidelberg: Springer Berlin Heidelberg, 1995, pp. 175–189.
  * [10] S. Rüping, “Learning with local models,” in _Local Pattern Detection_ , K. Morik, J.-F. Boulicaut, and A. Siebes, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2005, pp. 153–170.
  * [11] C. Mavridis and J. Baras, “Annealing optimization for progressive learning with stochastic approximation,” _arXiv preprint arXiv:2209.02826_ , 2022\.
  * [12] C. N. Mavridis, G. P. Kontoudis, and J. S. Baras, “Sparse gaussian process regression using progressively growing learning representations,” in _2022 61st IEEE Conference on Decision and Control (CDC). IEEE_ , 2022.
  * [13] C. N. Mavridis and J. S. Baras, “Progressive graph partitioning based on information diffusion,” in _IEEE Conference on Decision and Control_ , 2021, pp. 37–42.
  * [14] ——, “Online deterministic annealing for classification and clustering,” _IEEE Transactions on Neural Networks and Learning Systems_ , 2022.
  * [15] K. Rose, “Deterministic annealing for clustering, compression, classification, regression, and related optimization problems,” _Proceedings of the IEEE_ , vol. 86, no. 11, pp. 2210–2239, 1998.
  * [16] V. S. Borkar, _Stochastic approximation: a dynamical systems viewpoint_. Springer, 2009, vol. 48.
  * [17] A. Banerjee, S. Merugu, I. S. Dhillon, and J. Ghosh, “Clustering with bregman divergences,” _Journal of machine learning research_ , vol. 6, no. Oct, pp. 1705–1749, 2005.
  * [18] T. Villmann, S. Haase, F.-M. Schleif, B. Hammer, and M. Biehl, “The mathematics of divergence based online learning in vector quantization,” in _IAPR Workshop on Artificial Neural Networks in Pattern Recognition_. Springer, 2010, pp. 108–119.
  * [19] J. Bruna and S. Mallat, “Invariant scattering convolution networks,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 35, no. 8, pp. 1872–1886, 2013.
  * [20] C. N. Mavridis and J. S. Baras, “Convergence of stochastic vector quantization and learning vector quantization with bregman divergences,” _IFAC-PapersOnLine_ , vol. 53, no. 2, 2020.
  * [21] E. T. Jaynes, “Information theory and statistical mechanics,” _Physical review_ , vol. 106, no. 4, p. 620, 1957.
  * [22] C. Mavridis, E. Noorani, and J. S. Baras, “Risk sensitivity and entropy regularization in prototype-based learning,” in _2022 30th Mediterranean Conference on Control and Automation (MED)_. IEEE, 2022, pp. 194–199.
  * [23] L. Devroye, L. Györfi, and G. Lugosi, _A probabilistic theory of pattern recognition_. Springer Science & Business Media, 2013, vol. 31.
  * [24] R. M. Gray, “Vector quantization,” _Readings in speech recognition_ , vol. 1, no. 2, pp. 75–100, 1990.
  * [25] E. A. Riskin and R. M. Gray, “A greedy tree growing algorithm for the design of variable rate vector quantizers (image compression),” _IEEE Transactions on Signal Processing_ , vol. 39, no. 11, pp. 2500–2507, 1991.
  * [26] A. B. Nobel and R. A. Olshen, “Termination and continuity of greedy growing for tree-structured vector quantizers,” _IEEE Transactions on Information Theory_ , vol. 42, no. 1, pp. 191–205, 1996.
  * [27] S. Mallat, “Understanding deep convolutional networks,” _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , vol. 374, no. 2065, p. 20150203, 2016.
  * [28] ——, _A wavelet tour of signal processing_. Elsevier, 1999.
  * [29] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” _SIAM review_ , vol. 43, no. 1, pp. 129–159, 2001.
  * [30] Y. LeCun, “The next frontier in ai: Unsupervised learning,” _https://www.youtube.com/watch?v=IbjF5VjniVE_ , 2016.
  * [31] J. S. Baras and S. I. Wolk, “Efficient organization of large ship radar databases using wavelets and structured vector quantization,” in _Proceedings of 27th Asilomar Conference on Signals, Systems and Computers_. IEEE, 1993, pp. 491–498.
  * [32] ——, “Wavelet-based progressive classification of high-range resolution radar returns,” in _Wavelet Applications_ , H. H. Szu, Ed., vol. 2242, International Society for Optics and Photonics. SPIE, 1994, pp. 967 – 977. [Online]. Available: https://doi.org/10.1117/12.170034
  * [33] M. Andreux, T. Angles, G. Exarchakis, R. Leonarduzzi, G. Rochette, L. Thiry, J. Zarka, S. Mallat, J. Andén, E. Belilovsky, J. Bruna, V. Lostanlen, M. J. Hirn, E. Oyallon, S. Zhang, C. Cella, and M. Eickenberg, “Kymatio: Scattering transforms in python,” 2019.
  * [34] F. Anselmi, L. Rosasco, C. Tan, and T. Poggio, “Deep convolutional networks are hierarchical kernel machines,” _arXiv preprint arXiv:1508.01084_ , 2015\.
  * [35] S. Mallat, “Group invariant scattering,” _Communications on Pure and Applied Mathematics_ , vol. 65, no. 10, pp. 1331–1398, 2012.
  * [36] S. Milani, N. Topin, M. Veloso, and F. Fang, “A survey of explainable reinforcement learning,” _arXiv preprint arXiv:2202.08434_ , 2022.
  * [37] L. Breiman, “Random forests,” _Machine learning_ , vol. 45, no. 1, pp. 5–32, 2001.
  * [38] V. S. Borkar, “Stochastic approximation with two time scales,” _Systems & Control Letters_, vol. 29, no. 5, pp. 291–294, 1997.

##  Appendix A Proof of Lemma 6 (Derivation of the Association Probabilities).

Recall that by definition (4), we get

$\displaystyle F(\mu)$ $\displaystyle:=(1-\lambda)\int
p(x)\sum_{i}p(\mu_{i}|x)d(x,\mu_{i})~{}dx$ $\displaystyle\quad+\lambda\int
p(x)\sum_{i}p(\mu_{i}|x)\log p(\mu_{i}|x)~{}dx-\lambda H(X)$

We form the Lagrangian:

$\displaystyle\mathcal{L}_{f}$
$\displaystyle(\left\\{p(\mu_{i}|x)\right\\},\nu):=$ (39)
$\displaystyle=(1-\lambda)D(\mu)-\lambda
H(\mu)+\nu\left(\sum_{i}p(\mu_{i}|x)-1\right)$ $\displaystyle=(1-\lambda)\int
p(x)\sum_{i}p(\mu_{i}|x)d(x,\mu_{i})~{}dx$ $\displaystyle\quad+\lambda\int
p(x)\sum_{i}p(\mu_{i}|x)\log p(\mu_{i}|x)~{}dx$
$\displaystyle\quad+\nu\left(\sum_{i}p(\mu_{i}|x)-1\right)-\lambda\mathbb{E}\left[-\log
p(X)\right]$

Taking $\frac{\partial\mathcal{L}}{\partial p(\mu|x)}=0$ yields:

$\displaystyle(1-\lambda)d(x,\mu_{i})+\lambda(1+\log p(\mu_{i}|x))+\nu=0$
$\displaystyle\implies\log
p(\mu_{i}|x)=-\frac{1-\lambda}{\lambda}d(x,\mu_{i})-\left(1+\frac{\nu}{\lambda}\right)$
$\displaystyle\implies
p(\mu_{i}|x)=\frac{e^{-\frac{1-\lambda}{\lambda}d(x,\mu_{i})}}{e^{1+\frac{\nu}{\lambda}}}$

Finally, from the condition $\sum_{i}p(\mu_{i}|x)=1$, it follows that

$e^{1+\frac{\nu}{\lambda}}=\sum_{i}e^{-\frac{1-\lambda}{\lambda}d(x,\mu_{i})}$

which completes the proof.

##  Appendix B Proof of Theorem 3 (Convergence of the Online Learning Rule).

We are going to use fundamental results from stochastic approximation theory.
For completeness, we present the key theorems in what follows.

######  Theorem 10 ([16], Ch.2).

Almost surely, the sequence $\left\\{x_{n}\right\\}\in
S\subseteq\mathbb{R}^{d}$ generated by the following stochastic approximation
scheme:

$\displaystyle x_{n+1}=x_{n}+\alpha(n)\left[h(x_{n})+M_{n+1}\right],\ n\geq 0$
(40)

with prescribed $x_{0}$, converges to a (possibly sample path dependent)
compact, connected, internally chain transitive, invariant set of the o.d.e:

$\displaystyle\dot{x}(t)=h\left(x(t)\right),~{}t\geq 0,$ (41)

where $x:\mathbb{R}_{+}\rightarrow\mathbb{R}_{d}$ and $x(0)=x_{0}$, provided
the following assumptions hold:

  * (A1)

The map $h:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is Lipschitz in $S$, i.e.,
$\exists L$ with $0<L<\infty$ such that $\left\|h(x)-h(y)\right\|\leq
L\left\|x-y\right\|,~{}x,y\in S$,

  * (A2)

The stepsizes $\left\\{\alpha(n)\in\mathbb{R}_{++},~{}n\geq 0\right\\}$
satisfy $\sum_{n}\alpha(n)=\infty$, and $\sum_{n}\alpha^{2}(n)<\infty$ ,

  * (A3)

$\left\\{M_{n}\right\\}$ is a martingale difference sequence with respect to
the increasing family of $\sigma$-fields
$\mathcal{F}_{n}:=\sigma\left(x_{m},M_{m},~{}m\leq n\right)$, ${n\geq 0}$,
i.e., $\mathbb{E}\left[M_{n+1}|\mathcal{F}_{n}\right]=0~{}a.s.$, for all
$n\geq 0$, and $\left\\{M_{n}\right\\}$ are square-integrable with
$\mathbb{E}\left[\left\|M_{n+1}\right\|^{2}|\mathcal{F}_{n}\right]\leq
K\left(1+\left\|x_{n}\right\|^{2}\right),~{}a.s.$, where $n\geq 0$ for some
$K>0$,

  * (A4)

The iterates $\left\\{x_{n}\right\\}$ remain bounded a.s., i.e.,
${\sup_{n}\left\|x_{n}\right\|<\infty}$ $a.s.$

As an immediate result, the following corollary also holds:

######  Corollary 10.1.

If the only internally chain transitive invariant sets for (41) are isolated
equilibrium points, then, almost surely, $\left\\{x_{n}\right\\}$ converges to
a, possibly sample dependent, equilibrium point of (41).

Now we are in place to prove the following theorem:

######  Theorem 11.

Let $S$ a vector space, $\mu\in S$, and $X:\Omega\rightarrow S$ be a random
variable defined in a probability space
$\left(\Omega,\mathcal{F},\mathbb{P}\right)$. Let $\left\\{x_{n}\right\\}$ be
a sequence of independent realizations of $X$, and
$\left\\{\alpha(n)>0\right\\}$ a sequence of stepsizes such that
$\sum_{n}\alpha(n)=\infty$, and $\sum_{n}\alpha^{2}(n)<\infty$. Then the
random variable $m_{n}=\nicefrac{{\sigma_{n}}}{{\rho_{n}}}$, where
$(\rho_{n},\sigma_{n})$ are sequences defined by

$\displaystyle\rho_{n+1}$
$\displaystyle=\rho_{n}+\alpha(n)\left[p(\mu|x_{n})-\rho_{n}\right]$ (42)
$\displaystyle\sigma_{n+1}$
$\displaystyle=\sigma_{n}+\alpha(n)\left[x_{n}p(\mu|x_{n})-\sigma_{n}\right],$

converges to $\mathbb{E}\left[X|\mu\right]$ almost surely, i.e.
$m_{n}\xrightarrow{a.s.}\mathbb{E}\left[X|\mu\right]$.

###### Proof.

We will use the facts that $p(\mu)=\mathbb{E}\left[p(\mu|x)\right]$ and
$\mathbb{E}\left[\mathds{1}_{\left[\mu\right]}X\right]=\mathbb{E}\left[xp(\mu|x)\right]$.
The recursive equations (42) are stochastic approximation algorithms of the
form:

$\displaystyle\rho_{n+1}$
$\displaystyle=\rho_{n}+\alpha(n)[(p(\mu)-\rho_{n})+(p(\mu|x_{n})-\mathbb{E}\left[p(\mu|X)\right])]$
$\displaystyle\sigma_{n+1}$
$\displaystyle=\sigma_{n}+\alpha(n)[(\mathbb{E}\left[\mathds{1}_{\left[\mu\right]}X\right]-\sigma_{n})+$
$\displaystyle\quad\quad\quad\quad\quad\quad(x_{n}p(\mu|x_{n})-\mathbb{E}\left[x_{n}p(\mu|X)\right])]$

It is obvious that both stochastic approximation algorithms satisfy the
conditions of Theorem 10 and Corollary 10.1. As a result, they converge to the
asymptotic solution of the differential equations

$\displaystyle\dot{\rho}$ $\displaystyle=p(\mu)-\rho$
$\displaystyle\dot{\sigma}$
$\displaystyle=\mathbb{E}\left[\mathds{1}_{\left[\mu\right]}X\right]-\sigma$

which can be trivially derived through standard ODE analysis to be
$\left(p(\mu),\mathbb{E}\left[\mathds{1}_{\left[\mu\right]}X\right]\right)$.
In other words, we have shown that

$\left(\rho_{n},\sigma_{n}\right)\xrightarrow{a.s.}\left(p(\mu),\mathbb{E}\left[\mathds{1}_{\left[\mu\right]}X\right]\right)$

The convergence of $m_{n}$ follows from the fact that
$\mathbb{E}\left[X|\mu\right]=\nicefrac{{\mathbb{E}\left[\mathds{1}_{\left[\mu\right]}X\right]}}{{p(\mu)}}$,
and standard results on the convergence of the product of two random
variables. ∎

##  Appendix C Proof of Theorem 4 (Consistency of ODA as a Density Estimator).

According to Theorem 3, as $n\rightarrow\infty$, the stochastic approximation
algorithm in (12), (13) minimizes the cost function $F^{*}$ in (9). Moreover,
it is easy to see that, in the limit $\lambda\rightarrow 0$, we get

$\lim_{\lambda\rightarrow 0}p^{*}(\mu_{i}|x)=\lim_{\lambda\rightarrow
0}\frac{e^{-\frac{1-\lambda}{\lambda}d(x,\mu_{i})}}{\sum_{j}e^{-\frac{1-\lambda}{\lambda}d(x,\mu_{j})}}=\mathds{1}_{\left[x\in
S_{i}\right]}$

and

$\displaystyle\lim_{\lambda\rightarrow 0}F^{*}(\mu)=J(\mu)$
$\displaystyle=\mathbb{E}\left[\min_{i}d(X,\mu_{i})\right]$
$\displaystyle=\int p(x)\sum_{i}\mathds{1}_{\left[x\in
S_{i}\right]}d_{\phi}(x,\mu_{i})~{}dx$
$\displaystyle=\sum_{i}\int_{S_{i}}p(x)d_{\phi}(x,\mu_{i})~{}dx$

where $S_{i}=\left\\{x\in
S:i=\operatorname*{arg\,min}\limits_{j}~{}d(x,\mu_{j})\right\\}$. In addition,
due to the bifurcation phenomenon, $\lambda\rightarrow 0$, induces
$k\rightarrow\infty$.

Next we show that as the number of prototypes goes to infinity, i.e., if
$k\rightarrow\infty$, we get $\min_{\mu}J(\mu)=0$. First, consider a sub-
optimal solution $w:=\left\\{w_{j}\right\\}_{j=1}^{k}$, with
$p(w_{i}|x)=\mathds{1}_{\left[x\in\Sigma_{i}\right]}$ with the property that
$Vol(\Sigma_{i})=\int_{\Sigma_{i}}dx=O(\frac{1}{k})$, i.e., the Voronoi cells
$\Sigma_{i}(k)$ form a roughly uniform partition. In that case

$\lim_{k\rightarrow\infty}J(w)=\lim_{k\rightarrow\infty}\sum_{i}\int_{\Sigma_{i}}p(x)d_{\phi}(x,w_{i})~{}dx=0$

where we have used the continuity of the density, the compactness of $S$, and
the fact that $\lim_{k\rightarrow\infty}Vol(\Sigma_{j})=0$ since
$Vol(\Sigma_{j})=O(\frac{1}{k})$. We note that these convergence results hold
as long as $\frac{k}{n}\rightarrow 0$, i.e., the rate of increase of $k$ is
lower than that of the number of observations $n$.

As a result, since $0\leq J(\mu)\leq J(w)\rightarrow 0$, due to the optimality
of $\mu$, it follows that $J(\mu)\rightarrow 0$ a.s., as well. This implies
that $\lim_{k\rightarrow\infty}Vol(S_{i})\rightarrow 0$. Now define the random
variable

$Y_{n}:=\frac{\mathds{1}_{\left[x_{n}\in S_{i}\right]}}{Vol(S_{i})}$

where

$\mathbb{E}\left[Y_{n}\right]:=\frac{\mathbb{E}\left[\mathds{1}_{\left[x_{n}\in
S_{i}\right]}\right]}{Vol(S_{i})}=\frac{\mathbb{P}\left[x_{n}\in
S_{i}\right]}{Vol(S_{i})}=\frac{\int_{x\in S_{i}}p(x)~{}dx}{\int_{x\in
S_{i}}~{}dx}=\bar{Y}$ (43)

From the Strong Law of Large Numbers (SLLN), we get that

$\frac{1}{n}\frac{\sum_{n}\mathds{1}_{\left[x_{n}\in
S_{i}\right]}}{Vol(S_{i})}\rightarrow\bar{Y},\ a.s.$

The claim that for $n\rightarrow\infty$, $k\rightarrow\infty$, and
$\frac{k}{n}\rightarrow 0$,

$\hat{p}(x):=\frac{1}{n}\frac{\sum_{n}\mathds{1}_{\left[x_{n}\in
S_{i}\right]}}{Vol(S_{i})}\rightarrow p(x),\ a.s.$

follows by observing that $\lim_{k\rightarrow\infty}\bar{Y}=p(x)$ which
follows from (43) and the fact that the limit
$\lim_{k\rightarrow\infty}Vol(S_{i})\rightarrow 0,\ \forall i$.

##  Appendix D Proof of Theorem 5 (Proof of Bifurcation Phenomena).

To obtain the optimality condition (16) as a function of the temperature level
$\lambda$, we recall that

$\displaystyle F^{*}(y)=(1-\lambda)\int
p(x)\sum_{i}p(y_{i}|x)d_{\phi}(x,y_{i})~{}dx$ $\displaystyle\quad+\lambda\int
p(x)\sum_{i}p(y_{i}|x)\log p(y_{i}|x)~{}dx+\lambda\int p(x)\log p(X)$

where $y=\mu+\epsilon\psi$, and $d_{\phi}$ is a Bregman divergence for an
appropriately defined strictly convex function $\phi$. By direct
differentiation, we can compute

$\frac{d}{d\epsilon}d_{\phi}(x,y_{i})=-\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}(x-y_{i})^{\mathrm{T}}\psi$ (44)

and

$\frac{d^{2}}{d\epsilon^{2}}d_{\phi}(x,y_{i})=\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}\psi^{\mathrm{T}}\psi$ (45)

and given (6), we get

$\displaystyle\frac{d}{d\epsilon}$ $\displaystyle
p(y_{i}|x)=\frac{1-\lambda}{\lambda}\sum_{j}p(y_{i}|x)p(y_{j}|x)\bigg{[}\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}(x-y_{i})^{\mathrm{T}}\psi_{i}-$ (46)
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\frac{\partial^{2}\phi(y_{j})}{\partial
y_{j}^{2}}(x-y_{j})^{\mathrm{T}}\psi_{j}\bigg{]}$
$\displaystyle=\frac{1-\lambda}{\lambda}p(y_{i}|x)\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}(x-y_{i})^{\mathrm{T}}\psi_{i}$
$\displaystyle-\frac{1-\lambda}{\lambda}p(y_{i}|x)\sum_{j}p(y_{j}|x)\frac{\partial^{2}\phi(y_{j})}{\partial
y_{j}^{2}}(x-y_{j})^{\mathrm{T}}\psi_{j}$
$\displaystyle=-\frac{1-\lambda}{\lambda}p(y_{i}|x)\frac{d}{d\epsilon}d_{\phi}(x,y_{i})$
$\displaystyle+\frac{1-\lambda}{\lambda}p(y_{i}|x)\sum_{j}p(y_{j}|x)\frac{d}{d\epsilon}d_{\phi}(x,y_{j})$

Now the optimality condition takes the form:

$\displaystyle\frac{d^{2}}{d\epsilon^{2}}F^{*}(y)$
$\displaystyle=(1-\lambda)\int
p(x)\sum_{i}\frac{d^{2}}{d\epsilon^{2}}\left(p(y_{i}|x)d_{\phi}(x,y_{i})\right)~{}dx$
(47) $\displaystyle+\lambda\int
p(x)\sum_{i}\frac{d^{2}}{d\epsilon^{2}}\left(p(y_{i}|x)\log
p(y_{i}|x)\right)~{}dx$

Equation (47) uses the terms given in (48), and (49) below:

$\displaystyle\frac{d^{2}}{d\epsilon^{2}}$
$\displaystyle\left(p(y_{i}|x)d_{\phi}(x,y_{i})\right)=p(y_{i}|x)\frac{d^{2}}{d\epsilon^{2}}d_{\phi}(x,y_{i})$
(48)
$\displaystyle\quad+2\frac{d}{d\epsilon}d_{\phi}(x,y_{i})\frac{d}{d\epsilon}p(y_{i}|x)+d_{\phi}(x,y_{i})\frac{d^{2}}{d\epsilon^{2}}p(y_{i}|x)$
$\displaystyle\frac{d^{2}}{d\epsilon^{2}}$ $\displaystyle\left(p(y_{i}|x)\log
p(y_{i}|x)\right)=\frac{d^{2}}{d\epsilon^{2}}p(y_{i}|x)$ (49)
$\displaystyle+\log
p(y_{i}|x)\frac{d^{2}}{d\epsilon^{2}}p(y_{i}|x)+\frac{1}{p(y_{i}|x)}\left(\frac{d}{d\epsilon}p(y_{i}|x)\right)^{2}$

First, notice that
$\sum_{i}\frac{d^{2}}{d\epsilon^{2}}p(y_{i}|x)=\frac{d^{2}}{d\epsilon^{2}}\sum_{i}p(y_{i}|x)=0$.
Using the expressions (44), (45), and (46) we get (50), (51), (52), (53) that
read as:

$\sum_{i}p(y_{i}|x)\frac{d^{2}}{d\epsilon^{2}}d_{\phi}(x,y_{i})=\sum_{i}p(y_{i}|x)\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}\psi_{i}^{\mathrm{T}}\psi_{i}$ (50) $\displaystyle\sum_{i}2$
$\displaystyle\frac{d}{d\epsilon}d_{\phi}(x,y_{i})\frac{d}{d\epsilon}p(y_{i}|x)=$
(51)
$\displaystyle=-2\frac{1-\lambda}{\lambda}\sum_{i}p(y_{i}|x)\left(\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}(x-y_{i})^{\mathrm{T}}\psi_{i}\right)^{2}$
$\displaystyle+2\frac{1-\lambda}{\lambda}\sum_{i}p(y_{i}|x)\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}(x-y_{i})^{\mathrm{T}}\psi_{i}$
$\displaystyle\quad\quad\quad\quad\quad\sum_{j}p(y_{j}|x)\frac{\partial^{2}\phi(y_{j})}{\partial
y_{j}^{2}}(x-y_{j})^{\mathrm{T}}\psi_{j}$
$\displaystyle=\sum_{i}p(y_{i}|x)\left(\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}\right)^{2}$
$\displaystyle\quad\quad\psi_{i}^{\mathrm{T}}\left[-2\frac{1-\lambda}{\lambda}(x-y_{i})(x-y_{i})^{\mathrm{T}}\right]\psi_{i}$
$\displaystyle\quad+2\frac{1-\lambda}{\lambda}\left(\sum_{i}p(y_{i}|x)\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}(x-y_{i})^{\mathrm{T}}\psi_{i}\right)^{2}$ $\displaystyle\sum_{i}$
$\displaystyle\log p(y_{i}|x)\frac{d^{2}}{d\epsilon^{2}}p(y_{i}|x)=$ (52)
$\displaystyle=-\sum_{i}\frac{1-\lambda}{\lambda}d_{\phi}(x,y_{i})\frac{d^{2}}{d\epsilon^{2}}p(y_{i}|x)$
$\displaystyle-\sum_{i}\log\left(\sum_{j}e^{-\frac{1-\lambda}{\lambda}d(x,y_{j})}\right)\frac{d^{2}}{d\epsilon^{2}}p(y_{i}|x)$
$\displaystyle=-\sum_{i}\frac{1-\lambda}{\lambda}d_{\phi}(x,y_{i})\frac{d^{2}}{d\epsilon^{2}}p(y_{i}|x)$
$\displaystyle-\left(\log\sum_{j}e^{-\frac{1-\lambda}{\lambda}d(x,y_{j})}\right)\sum_{i}\frac{d^{2}}{d\epsilon^{2}}p(y_{i}|x)$
$\displaystyle=-\frac{1-\lambda}{\lambda}\sum_{i}d_{\phi}(x,y_{i})\frac{d^{2}}{d\epsilon^{2}}p(y_{i}|x)$
$\displaystyle\sum_{i}$
$\displaystyle\frac{1}{p(y_{i}|x)}\left(\frac{d}{d\epsilon}p(y_{i}|x)\right)^{2}=$
(53)
$\displaystyle=\frac{\left(1-\lambda\right)^{2}}{\lambda^{2}}\sum_{i}p(y_{i}|x)\left(\frac{d}{d\epsilon}d_{\phi}(x,y_{i})\right)^{2}$
$\displaystyle+\frac{\left(1-\lambda\right)^{2}}{\lambda^{2}}\sum_{i}p(y_{i}|x)\left(\sum_{y_{j}}p(y_{j}|x)\frac{d}{d\epsilon}d_{\phi}(x,y_{j})\right)^{2}$
$\displaystyle-2\frac{\left(1-\lambda\right)^{2}}{\lambda^{2}}\left(\sum_{i}p(y_{i}|x)\frac{d}{d\epsilon}d_{\phi}(x,y_{i})\right)^{2}$
$\displaystyle=\frac{\left(1-\lambda\right)^{2}}{\lambda^{2}}\sum_{i}p(y_{i}|x)\left(\frac{d}{d\epsilon}d_{\phi}(x,y_{i})\right)^{2}$
$\displaystyle-\frac{\left(1-\lambda\right)^{2}}{\lambda^{2}}\left(\sum_{i}p(y_{i}|x)\frac{d}{d\epsilon}d_{\phi}(x,y_{i})\right)^{2}$

Finally, plugging (50), (51), (52), (53) in (47), the optimality condition
(16) becomes

$\displaystyle 0$ $\displaystyle=(1-\lambda)\int
p(x)\sum_{i}p(y_{i}|x)\frac{d^{2}}{d\epsilon^{2}}d_{\phi}(x,y_{i})dx$
$\displaystyle-\frac{\left(1-\lambda\right)^{2}}{\lambda}\int
p(x)\sum_{i}p(y_{i}|x)\left(\frac{d}{d\epsilon}d_{\phi}(x,y_{i})\right)^{2}dx$
$\displaystyle+\frac{\left(1-\lambda\right)^{2}}{\lambda}\int
p(x)\left(\sum_{i}p(y_{i}|x)\frac{d}{d\epsilon}d_{\phi}(x,y_{i})\right)^{2}dx$

which can be written more precisely as

$\displaystyle 0$ $\displaystyle=\int
p(x)\sum_{i}p(y_{i}|x)\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}\psi^{\mathrm{T}}$
$\displaystyle\left[I-\frac{1-\lambda}{\lambda}\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}(x-y_{i})(x-y_{i})^{\mathrm{T}}\right]\psi dx$
$\displaystyle+\frac{1-\lambda}{\lambda}\int
p(x)\left(\sum_{i}p(y_{i}|x)\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}(x-y_{i})^{\mathrm{T}}\psi\right)^{2}dx$

and finally as

$\displaystyle 0$
$\displaystyle=\sum_{i}p(y_{i})\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}\psi^{\mathrm{T}}\left[I-\frac{1-\lambda}{\lambda}\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}C_{x|y_{i}}\right]\psi$ (54)
$\displaystyle+\frac{1-\lambda}{\lambda}\int
p(x)\left(\sum_{i}p(y_{i}|x)\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}(x-y_{i})^{\mathrm{T}}\psi\right)^{2}dx$

where

$\displaystyle C_{x|y_{i}}:$
$\displaystyle=\mathbb{E}\left[(x-y_{i})(x-y_{i})^{\mathrm{T}}|y_{i}\right]$
$\displaystyle=\int p(x|y_{i})(x-y_{i})(x-y_{i})^{\mathrm{T}}dx$

The left-hand side of (54) is positive for all perturbations
$\left\\{\psi\right\\}$ if and only if the first term is positive. To see
that, notice that the second term of (54) is clearly non-negative. For the
left-hand side to be non-positive, the first term needs to be non-positive as
well, i.e., there should exist at least one codevector value, say $y_{n}$,
such that $p(y_{n})>0$ and
$\left[I-\frac{1-\lambda}{\lambda}\frac{\partial^{2}\phi(y_{n})}{\partial
y_{n}^{2}}C_{x|y_{n}}\right]\preceq 0$. In this case, there always exist a
perturbation vector $\left\\{y\right\\}$ such that $y=0$, $\forall y\neq
y_{n}$, and $\sum_{y=y_{n}}\psi=0$, that vanishes the second term, i.e.,
$\frac{1-\lambda}{\lambda}\int
p(x)\left(\sum_{i}p(y_{i}|x)\frac{\partial^{2}\phi(y_{i})}{\partial
y_{i}^{2}}(x-y_{i})^{\mathrm{T}}\psi\right)^{2}dx$ $=0$. In other words we
have shown that

$\displaystyle\frac{d^{2}}{d\epsilon^{2}}F^{*}(y)$
$\displaystyle>0\quad\Leftrightarrow\quad\exists y_{n}\text{ s.t. }p(y_{n})>0$
$\displaystyle\text{ and
}\left[I-\frac{1-\lambda}{\lambda}\frac{\partial^{2}\phi(y_{n})}{\partial
y_{n}^{2}}C_{x|y_{n}}\right]\succ 0$

which means that bifurcation occurs under the following condition

$\exists y_{n}\text{ s.t. }p(y_{n})>0\text{ and
}\det\left[I-\frac{1-\lambda}{\lambda}\frac{\partial^{2}\phi(y_{n})}{\partial
y_{n}^{2}}C_{x|y_{n}}\right]=0$

which completes the proof.

##  Appendix E Proof of Theorem 6 (Convergence of the Two-Timescale Online
Learning Rule).

It follows directly from the following fundamental result from stochastic
approximation theory that is included for completeness:

######  Theorem 12 (Ch. 6 of [38]).

Consider the sequences $\left\\{x_{n}\right\\}\in S\subseteq\mathbb{R}^{d}$
and $\left\\{y_{n}\right\\}\in\Sigma\subseteq\mathbb{R}^{k}$, generated by the
iterative stochastic approximation schemes:

$\displaystyle
x_{n+1}=x_{n}+\beta(n)\left[f(x_{n},y_{n})+M_{n+1}^{(x)}\right]$ (55)
$\displaystyle
y_{n+1}=y_{n}+\alpha(n)\left[g(x_{n},y_{n})+M_{n+1}^{(y)}\right]$ (56)

for $n\geq 0$ and $M_{n}^{(x)}$, $M_{n}^{(y)}$ martingale difference
sequences, and assume that $\sum_{n}\alpha(n)=\sum_{n}\beta(n)=\infty$,
$\sum_{n}(\alpha^{2}(n)+\beta^{2}(n))<\infty$, and
$\nicefrac{{\alpha(n)}}{{\beta(n)}}\rightarrow 0$, with the last condition
implying that the iterations for $\left\\{y_{n}\right\\}$ run on a slower
timescale than those for $\left\\{x_{n}\right\\}$. If the equation

$\dot{x}(t)=f(x(t),y),\ x(0)=x_{0}$

has an asymptotically stable equilibrium $\lambda(y)$ for fixed $y$ and some
Lipschitz mapping $\lambda$, and the equation

$\dot{y}(t)=g(\lambda(y(t)),y(t)),\ y(0)=y_{0}$

has an asymptotically stable equilibrium $y^{*}$, then, almost surely,
$(x_{n},y_{n})$ converges to $(\lambda(y^{*}),y^{*})$.

##  Appendix F Proof of Theorem 7 (Convergence to Bayes-Optimal Classifier).

The Bayes risk for minimum probability of error is given by

$\displaystyle J_{B}(\mu,c_{\mu}):$
$\displaystyle=\pi_{1}\sum_{i:c_{\mu_{i}}=0}\mathbb{P}\left\\{X\in
S_{i}|c=1\right\\}$ (57)
$\displaystyle+\pi_{0}\sum_{i:c_{\mu_{i}}=1}\mathbb{P}\left\\{X\in
S_{i}|c=0\right\\}$ (58)

where $\mathbb{P}\left\\{X\in S_{i}|c=j\right\\}=\int_{S_{i}}p(x|c=j)~{}dx$,
and $p(x|c=j)$ represents the conditional probability density function.

Consider the observations $\left\\{(x_{n},c_{n})\right\\}$, and let
$\hat{p}(x|c=j)$, $j=1,2$, be strongly consistent density estimators. Then the
estimated risk

$\displaystyle\hat{J}_{B}(\mu,c_{\mu}):$
$\displaystyle=\hat{\pi}_{1}\sum_{i:c_{\mu_{i}}=0}\mathbb{\hat{P}}\left\\{X\in
S_{i}|c=1\right\\}$ (59)
$\displaystyle+\hat{\pi}_{0}\sum_{i:c_{\mu_{i}}=1}\mathbb{\hat{P}}\left\\{X\in
S_{i}|c=0\right\\}$ (60)

where $\mathbb{\hat{P}}\left\\{X\in
S_{i}|c=j\right\\}=\int_{S_{i}}\hat{p}(x|c=j)$ and
$\hat{\pi}_{j}=\frac{\sum_{n}\mathds{1}_{\left[c_{n}=j\right]}}{n}$, converges
almost surely to $J_{B}$, i.e.,

$\hat{J}_{B}(\mu,c_{\mu})\rightarrow J_{B}(\mu,c_{\mu}),\ a.s.$ (61)

as $n\rightarrow\infty$. This follows from the fact that

$\hat{\pi}_{j}\hat{p}(x|c=j)\rightarrow\pi_{j}p(x|c=j),\ a.s.$ (62)

and the Lebesgue dominated convergence theorem. Therefore, the classification
rule $\hat{c}=\operatorname*{arg\,max}_{j}\hat{\pi}_{j}\hat{p}(x|c=j)$
converges to a Bayes-optimal classification rule.

|  Christos N. Mavridis (M’20) received the Diploma degree in electrical and
computer engineering from the National Technical University of Athens, Greece,
in 2017, and the M.S. and Ph.D. degrees in electrical and computer engineering
at the University of Maryland, College Park, MD, USA, in 2021. His research
interests include learning theory, stochastic optimization, systems and
control theory, multi-agent systems, and robotics. He is currently a
postdoctoral associate at the University of Maryland, and a visiting
postdoctoral fellow at KTH Royal Institute of Technology, Stockholm. He has
worked as a research intern for the Math and Algorithms Research Group at
Nokia Bell Labs, NJ, USA, and the System Sciences Lab at Xerox Palo Alto
Research Center (PARC), CA, USA. Dr. Mavridis is an IEEE member, and a member
of the Institute for Systems Research (ISR) and the Autonomy, Robotics and
Cognition (ARC) Lab. He received the Ann G. Wylie Dissertation Fellowship in
2021, and the A. James Clark School of Engineering Distinguished Graduate
Fellowship, Outstanding Graduate Research Assistant Award, and Future Faculty
Fellowship, in 2017, 2020, and 2021, respectively. He has been a finalist in
the Qualcomm Innovation Fellowship US, San Diego, CA, 2018, and he has
received the Best Student Paper Award (1st place) in the IEEE International
Conference on Intelligent Transportation Systems (ITSC), 2021.
---|---
|  John S. Baras (LF’13) received the Diploma degree in electrical and
mechanical engineering from the National Technical University of Athens,
Athens, Greece, in 1970, and the M.S. and Ph.D. degrees in applied mathematics
from Harvard University, Cambridge, MA, USA, in 1971 and 1973, respectively.
He is a Distinguished University Professor and holds the Lockheed Martin Chair
in Systems Engineering, with the Department of Electrical and Computer
Engineering and the Institute for Systems Research (ISR), at the University of
Maryland College Park. From 1985 to 1991, he was the Founding Director of the
ISR. Since 1992, he has been the Director of the Maryland Center for Hybrid
Networks (HYNET), which he co-founded. His research interests include systems
and control, optimization, communication networks, applied mathematics,
machine learning, artificial intelligence, signal processing, robotics,
computing systems, security, trust, systems biology, healthcare systems,
model-based systems engineering. Dr. Baras is a Fellow of IEEE (Life), SIAM,
AAAS, NAI, IFAC, AMS, AIAA, Member of the National Academy of Inventors and a
Foreign Member of the Royal Swedish Academy of Engineering Sciences. Major
honors include the 1980 George Axelby Award from the IEEE Control Systems
Society, the 2006 Leonard Abraham Prize from the IEEE Communications Society,
the 2017 IEEE Simon Ramo Medal, the 2017 AACC Richard E. Bellman Control
Heritage Award, the 2018 AIAA Aerospace Communications Award. In 2016 he was
inducted in the A. J. Clark School of Engineering Innovation Hall of Fame. In
2018 he was awarded a Doctorate Honoris Causa by his alma mater the National
Technical University of Athens, Greece.
---|---