Online Deterministic Annealing for Classification and Clustering

Unsupervised analysis can provide valuable insights into the nature of the dataset at hand, and it plays an important role in the context of visualization. Central to unsupervised learning is the representation of data in a vector space by typical representatives, which is formally defined in the following optimization problem:

When the training data are not available a priori but are being observed online, or when the processing of the entire dataset in every optimization iteration is computationally infeasible, a stochastic vector quantization algorithm can be defined as a recursive asynchronous stochastic approximation algorithm based on gradient descent [14]:

The supervised counterpart of vector quantization is the particularly attractive and intuitive approach of the competitive-learning Learning Vector Quantization (LVQ) algorithm, initially proposed by Kohonen [12]. LVQ for binary classification is formulated in the following optimization problem (and generalized to any type of classification task, see, e.g. [28]):

Prototype-based algorithms rely on measuring the proximity between different vector representations. In most cases the Euclidean distance or another convex metric is used, but this can be generalized to alternative dissimilarity measures inspired by information theory and statistical analysis, such as the Bregman divergences:

Notice that, as a divergence measure, Bregman divergence can be used to measure the dissimilarity of one probability distribution to another on a statistical manifold, and is a notion weaker than that of the distance. In particular, it does not need to be symmetric or satisfy the triangle inequality. In this work, we will concentrate on nonempty, compact convex sets $S\subseteq\mathbb{R}^{d}$ so that the derivative of $d_{\phi}$ with respect to the second argument can be written as

where $x,\mu\in S$, $\frac{\partial}{\partial\mu}$ represents differentiation with respect to the second argument of $d_{\phi}$, and $\nabla^{2}\phi(\mu)$ represents the Hessian matrix of $\phi$ at $\mu$.

The family of Bregman divergences provides proximity measures that have been shown to enhance the performance of a learning algorithm [31]. In addition, the following theorem shows that the use of Bregman divergences is both necessary and sufficient for the optimizer $\mu_{h}$ of (1) to be analytically computed as the expected value of the data inside $S_{h}$, which is implicitly used by many “centroid” algorithms, such as $k$-means [27]:

In the clustering problem (Problem 1), the distortion function $J$ is typically non convex and riddled with poor local minima. To partially deal with this phenomenon, soft-clustering approaches have been proposed as a probabilistic framework for clustering. In this case, an input vector $X$ is assigned, through the quantizer $Q$, to all codevectors $\mu_{h}\in M$ with probabilities $p(\mu_{h}|X)$ , where $\sum_{h=1}^{K}p(\mu_{h}|X)=1$. In this regard, the quantizer $Q:S\rightarrow M$ becomes a discrete random variable, with the set $M$ being its image, and can be fully described by the values of $M=\left\{\mu_{h}\right\}_{h=1}^{K}$ and the probability functions $\left\{p(\mu_{h}|x)\right\}_{h=1}^{K}$. In contrast, hard clustering assumes that $Q$ is a simple random variable that can be described fully by $M$ and $V=\left\{S_{h}\right\}_{h=1}^{K}$, since $p(\mu_{h}|X)=\mathds{1}_{\left[X\in S_{h}\right]}$ (see Problem 1).

For the randomized partition we can rewrite the expected distortion as

where $p(\mu|x)$ is the association probability relating the input vector $x$ with the codevector $\mu$. We note that, at the limit, where each input vector is assigned to a unique codevector with probability one, this reduces to the hard clustering distortion. The main idea in deterministic annealing, is to seek the distribution that minimizes $D$ subject to a specified level of randomness, measured by the Shannon entropy

by appealing to Jaynes’ maximum entropy principle¹¹1Informally, Jaynes’ principle states: of all the probability distributions that satisfy a given set of constraints, choose the one that maximizes the entropy. [32]. This multi-objective optimization is conveniently formulated as the minimization of the Lagrangian

where $T$ is the temperature parameter that acts as a Lagrange multiplier. Clearly, for large values of $T$ we maximize the entropy, and, as $T$ is lowered, we trade entropy for reduction in distortion. Equation (2) also represents the scalarization method for trade-off analysis between two performance metrics [33]. As $T$ varies we essentially transition from one Pareto point to another, and the sequence of the solutions will correspond to a Pareto curve of the multi-objective optimization (2) that resembles annealing processes in chemical engineering. In this regard, the entropy $H$, which is closely related to the “purity” of the clusters, acts as a regularization term which is given progressively less weight as $T$ decreases.

As in the case of vector quantization, we form a coordinate block optimization algorithm to minimize $F$, by successively minimizing it with respect to the association probabilities $p(\mu|x)$ and the codevector locations $\mu$. Minimizing $F$ with respect to the association probabilities $p(\mu|x)$ is straightforward and yields the Gibbs distribution

while, in order to minimize $F$ with respect to the codevector locations $\mu$ we set the gradients to zero

In the following theorem, we show that we can have analytical solution to the last optimization step (4) in a convenient centroid form, if $d$ is a Bregman divergence. This is a similar result to Theorem 1 for vector quantization.

This optimization procedure takes place for decreasing values of the temperature coefficient $T$ such that the solution maintains minimum free energy (thermal equilibrium) while gradually lowering the temperature. Adding to the physical analogy, it is significant that, as the temperature is lowered, the system undergoes a sequence of “phase transitions”, which consists of natural cluster splits where the cardinality of the codebook (number of clusters) increases. This is a bifurcation phenomenon and provides a useful tool for controlling the size of the clustering model relating it to the scale of the solution.

At very high temperature ($T\rightarrow\infty$) the optimization yields uniform association probabilities

and, provided $d:=d_{\phi}$ is a Bregman divergence, all the codevectors are located at the same point:

which is the expected value of $X$ (Theorem 1). This is true regardless of the number of codevectors available. We refer to the number of different codevectors resulting from the optimization process as effective codevectors. These define the cardinality of the codebook, which changes as we lower the temperature. The bifurcation occurs when the solution above a critical temperature $T_{c}$ is no longer the minimum of the free energy $F$ for $T<T_{c}$. A set of coincident codevectors then splits into separate subsets. These critical temperatures $T_{c}$ can be traced when the Hessian of $F$ loses its positive definite property, and are, in some cases, computable (see Theorem 1 in [17]). In other words, an algorithmic implementation needs only as many codevectors as the number of effective codevectors, which depends only on the temperature parameter, i.e. the Lagrange multiplier of the multi-objective minimization problem in (2). As will be shown in Section III-E, we can detect the bifurcation points by maintaining and perturbing pairs of codevectors at each effective cluster so that they separate only when a critical temperature is reached.

The conditional expectation $\mathbb{E}\left[X|\mu\right]$ in eq. (5) can be approximated by the sample mean of the data points weighted by their association probabilities $p(\mu|x)$, i.e.,

This approach, however, defines an offline (batch) optimization algorithm and requires the entire dataset to be available a priori, subtly assuming that it is possible to store and also quickly access the entire dataset at each iteration. This is rarely the case in practical applications and results to computationally costly iterations that are slow to converge. We propose an Online Deterministic Annealing (ODA) algorithm, that dynamically updates its estimate of the effective codevectors with every observation. This results in a significant reduction in complexity, that comes in two levels. The first refers to huge reduction in memory complexity, since we bypass the need to store the entire dataset, as well as the association probabilities $\left\{p(\mu|x),\ \forall x\right\}$ that map each data point in the dataset to each cluster. The second level refers to the nature of the optimization iterations. In the online approach the optimization iterations increase in number but become much faster, and practical convergence is often after a smaller number of observations.

To define an online training rule for the above optimization framework, we formulate a stochastic approximation algorithm to recursively estimate $\mathbb{E}\left[X|\mu\right]$ directly. Stochastic approximation, first introduced in [34], was originally conceived as a tool for statistical computation, and, since then, has become a central tool in a number of different disciplines, often times unbeknownst to the users, researchers and practitioners. It offers an online, adaptive, and computationally inexpensive optimization framework, properties that make it an ideal optimization method for machine learning algorithms. In addition to its connection with optimization and learning algorithms, however, stochastic approximation is strongly connected to dynamical systems, as well, a property that allows the study of its convergence through the analysis of an ordinary differential equation, as illustrated in the following theorem:

As an immediate result, the following corollary also holds:

Now we are in place to prove the following theorem:

As a direct consequence of this theorem, the following corollary provides an online learning rule that solves the optimization problem of the deterministic annealing algorithm.

Finally, the learning rule (12), (13) can be used to define a consistent (histogram) density estimator at the limit $T\rightarrow 0$. This follows from the fact that as $T\rightarrow 0$, the number of clusters $K$ goes to infinity, $p(\mu_{h}|X)\rightarrow\mathds{1}_{\left[X\in S_{h}\right]}$, and, as a result, $F\rightarrow J$, i.e., the consistency of Alg. 1 can be studied with similar arguments to the stochastic divergence-based vector quantization algorithm (1) (see [14, 13]).

We can extend the proposed learning algorithm to be used for classification as well. In this case we can rewrite the expected distortion as

where $d^{b}(c_{x},c_{\mu})=\mathds{1}_{\left[c_{x}\neq c_{\mu}\right]}$. Because $d^{b}$ is not differentiable, using similar principles as in the case of LVQ, we can instead approximate the optimal solution by solving the minimization problem for the following distortion measure

This particular choice for the distortion measure $d^{c}$ will lead to some interesting regularization properties of the proposed online approach (see Section III-E).

It is easy to show that the coordinate block optimization steps (3) and (5), in this case become:

respectively. In the last step, we have assumed that the class $c_{\mu}$ of each centroid $\mu$ is given and cannot be changed dynamically by the algorithm, which results to the minimization with respect to $\mu$ only. In a similar fashion, it can be shown that the online learning rule that solves the optimization problem of the deterministic annealing algorithm for classification, based on the distortion measure (14), is given by:

where

At the limit $T\rightarrow 0$, the quantization scheme described above equipped with a majority-vote classification rule is strongly Bayes risk consistent, i.e., converges to the optimal (Bayes) probability of error (see Ch. 21 in [13]). However, due to the choice of the distortion measure $d^{c}$ in (14) used in ODA for classification, the algorithm can be used to estimate consistent class-conditional density estimators, which define the natural classification rule:

where $h^{*}=\operatorname*{arg\,max}\limits_{\tau=1,\ldots,K}~{}p(\mu_{\tau}|x),~{}h\in\left\{1,\ldots,K\right\}$.

The proposed Online Deterministic Annealing (ODA) algorithm (Algorithm 1), is based on (15), (16), and can be used for both clustering and classification alike, depending on whether the data belong to a single (clustering) or several classes (classification).

Temperature Schedule. The temperature schedule $T_{i}$ plays an important role in the behavior of the algorithm. Starting at high temperature $T_{max}$ ensures the correct operation of the algorithm. The value of $T_{max}$ depends on the domain of the data and should be large enough such that there is only one effective codevector at $T=T_{max}$. When the range of the domain of the data is not known a priori, overestimation is recommended. The stopping temperature $T_{min}$ can be set a priori or be decided online depending on the performance of the model at each temperature level. The temperature step $dT_{i}=T_{i-1}-T_{i}$ should be small enough such that no critical temperature is missed. On the other hand, the smaller the step $dT_{i}$, the more optimization problems need to be solved. It is common practice to use the geometric series $T_{i+1}=\gamma T_{i}$.

Stochastic Approximation. Regarding the stochastic approximation stepsizes, simple time-based learning rates, e.g. of the form $\alpha_{n}=\nicefrac{{1}}{{a+bn}}$, have been sufficient for fast convergence in all our experiments so far. Convergence is checked with the condition $d_{\phi}(\mu_{i}^{n},\mu_{i}^{n-1})<\epsilon_{c}$ for a given threshold $\epsilon_{c}$ that can depend on the domain of $X$. Exploring adaptive learning rates would be an interesting research direction for the future.

Bifurcation and Perturbations. To every temperature level $T_{i}$, corresponds a set of effective codevectors $\left\{\mu_{j}\right\}_{j=1}^{K_{i}}$, which consist of the different solutions of the optimization problem (2) at $T_{i}$. Bifurcation, at $T_{i}$, is detected by maintaining a pair of perturbed codevectors $\left\{\mu_{j}+\delta,\mu_{j}-\delta\right\}$ for each effective codevector $\mu_{j}$ generated at $T_{i-1}$, i.e. for $j=1\ldots,K_{i-1}$. Using arguments from variational calculus [17], 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. In case of a merge, one of the perturbed codevectors is removed from the model. Therefore, the cardinality of the model is at most doubled at every temperature level. For classification, a perturbed codevector for each distinct class is generated.

Regularization. Merging is detected by the condition $d_{\phi}(\mu_{j},\mu_{i})<\epsilon_{n}$, where $\epsilon_{n}$ is a design parameter that acts as a regularization term for the model. Large values for $\epsilon_{n}$ (compared to the support of the data $X$) lead to fewer effective codevectors, while small $\epsilon_{n}$ values lead to a fast growth in the model size, which is connected to overfitting. It is observed that, for practical convergence, the perturbation noise $\delta$ is best to not exceed $\epsilon_{n}$. An additional regularization mechanism that comes as a natural consequence of the stochastic approximation learning rule, is the detection of idle codevectors. To see that, notice that the sequence $\rho_{i}(n)$ resembles an approximation of the probability $p(\mu_{i},c_{\mu_{i}})$. In the updates (12), (13), $\rho_{i}(n)$ becomes negligible ($\rho_{i}(n)<\epsilon_{r}$) if not updated by any nearby observed data, which is a natural criterion for removing the codevector $\mu_{i}$. This happens if all observed data samples $x_{n}$ are largely dissimilar to $\mu_{i}$. In classification, because of the choice of $d^{c}$ in (14), codevectors $\mu_{i}$ that are not assigned the same class as the data in their vicinity, will end up to be removed, as well. The threshold $\epsilon_{r}$ is a parameter that usually takes values near zero.

Complexity. The worst case complexity of Algorithm 1 behaves as $O(\sigma_{max}NK_{max}^{2}d),$ where:

• •

  $N$ is an upper bound of the number of data samples observed, which should be large enough to overestimate the iterations needed for convergence;

• •

  $d$ is the dimension of the input vectors, i.e., $x\in\mathbb{R}^{d}$;

• •

  $K_{max}$ is the maximum number of codevectors allowed;

• •

  $\sigma_{max}=\left\{\sigma_{1},\sigma_{2},\ldots,\sigma_{K_{max}}\right\}$, where $\sigma_{i}$ is the number of temperature values in our temperature schedule that lie between two critical temperatures $T_{i}$ and $T_{i+1}$, with the understanding that at $T_{i}$ there are $i$ distinct effective codevectors present. Here we have assumed that $K_{max}$ is achievable within our temperature schedule.

Fine-Tuning. In practice, because the convergence to the Bayes decision surface comes at the limit $(K,T)\rightarrow(\infty,0)$, a fine-tuning mechanism should be designed to run on top of the proposed algorithm after $T_{min}$. This can be either an LVQ algorithm (Section II-B) or some other local model.