Mean field models for large data–clustering problems

1.   Introduction

Particle and kinetic models for consensus and cluster formation appeared in recent literature for self–organized socio–economic dynamical systems as opinion formation, flocking of birds or fish, elections and referendums under influence of mass media, etc. See e.g. [5,11,18,19,20,23,25,34,16], the review articles [35,4,1] and the book [40]. A related research direction is based on using the consensus features of these models in an artificial way to solve problems of optimization or segmentation of data in large dimensions [28,41,15,29].

In this paper, we aim at formulating suitable models on the microscopic (or particle) level as well as on the mean–field (or kinetic) level to describe the partition of a large set of data in clusters. This problem is also known as data clustering problem and it is widely studied in many applications like pattern recognition, shape detection and image segmentation problems. The proposed methods do not need to have fixed a–priori number of clusters and clusters are characterized by small in–group and large out–group distances.

Since we are interested in modes characterizing distance and qualitative features without additional physical or socio–economical modeling background the proposed model will generalize the Hegselmann–Krause (HK) opinion dynamics model [25]. Originally, the HK model was proposed in a microscopic setting, in one–dimensional spatial and time discrete framework. Several extensions exist, in the sequel we will briefly review the basic model before discussing its extension towards the image clustering problems.

To this aim, let us consider a group of $n$ particles with a (scalar) initial state $x_i(0)\in \mathbb{R}$, $i = 1,\dots,n$, and the state of each particle varies depending on the state of the others. The key idea of the Hegselmann–Krause (HK) model [25] is that particles with completely different opinions do not influence each other, and a sort of mediation occurs among agents whose opinions are within a bounded confidence interval described by a parameter $\epsilon\geq0.$ Let ${\bf{{x}}}(t) = [x_1(t),\dots,x_n(t)]^T$ be the state of the system at time $t\geq 0$. Then the dynamic of the $i$–th particle is given by

where ${\bf{{A}}}(t,\epsilon)\in \mathbb{R}^{n\times n}$ is the time–varying adjacency matrix whose entries are in the form

with

defining the neighborhood of the $i$–th particle at time $t$, and

the type of interactions. Precisely, when $\sigma_i = n$ the interactions are symmetric since $A_{ij} = A_{ji}$, $\forall\,i,j$, but ${\bf{{A}}}$ is not a stochastic matrix. Instead, when $\sigma_i = \left|{ \mathcal{N}_i(t,\epsilon)}\right|$, then ${\bf{{A}}}$ is a right stochastic matrix but interactions are no longer symmetric.

Several works have been proposed in the literature which analyze the properties of the HK model. For instance, for the analysis in the time discrete setting we refer to [27]. In [7] it is proven that, during the evolution of the system (1), the order of the states is preserved. Thanks to the definition of the interaction kernel (2)–(4), if $\left|{x_i (t) - x_{i+1} (t)}\right| > \epsilon$, at some time $t$, it remains true for larger times. Therefore, the HK model tends to group the initial states in a finite number of clusters as proved in [35]. Following [35], we define a cluster $\mathcal{C}(t)$ at time $t\geq 0$ a subset of particles separated from all the other particles

In [7,22,30] the stability of the dynamical model is investigated. In particular, the fact that the system converges to a steady profile in finite time is proved in [7]. For further results on the one–dimensional local and symmetric model we refer also to [8]. We point out that also behavior of cluster formation in the transient is of interest in the mathematical literature [21].

In [36,10], the one-dimensional Hegselmann–Krause is generalized to the case of a multi–dimensional data–set. Subsequently, the multi–dimensional HK model has been used as a technique to cluster a big amount of data into a small number of subsets with some common features [38] and to compare its performance with the $k$–means algorithm [31]. Recently, in [29,37] new approaches to clustering problems and image segmentation have been proposed based on the Kuramoto model.

Here, we introduce a generalization of the multi–dimensional formulation of the HK model to solve data clustering problems. This amounts to take into account clustering with respect to different features. We deal with data having both time dependent and static features. The latter describe intrinsic properties of a datum, such as the measure of a trustworthy information or the color intensity of pixels in images. We derive the corresponding mean–field limit and investigate analytically the properties of the model. In particular, following [2] the mean–field formulation allows the use of a random subset algorithm for efficient computations of the clusters.

The rest of the manuscript is organized as follows. The microscopic model is introduced at the beginning of Section 2 and briefly discussed in Section 2.1. The proposed model is still a microscopic model and we describe the case of large data using a mean–field equation in Section 3. Analytical properties of the kinetic equation, such as a–priori estimation on the evolution of the moments and characterization of the limit distribution, are discussed in Section 3.1. Numerical evidence of the theoretical results is provided in Section 4.1. Further, we propose applications on detection and compression of data, such as shape detection, in Section 4.2, and image segmentation, in Section 4.3. We finally conclude with some remarks and future research directions in Section 5.

2.   Microscopic models for data–clustering

Each particle $i = 1,\dots,N$ is endowed with a time–dependent state vector ${\bf{{x}}}_i = {\bf{{x}}}_i(t)\in \mathbb{R}^{d_1}$ as well as features ${\bf{{c}}}_i = [c_{i,1},\dots,c_{i,d_2}] \in \mathbb{R}^{d_2}$ representing static characteristics of the system, i.e., ${\bf{{c}}}_i$ is independent of time. As a motivation example consider an image segmentation problem where ${\bf{{x}}}_i$ are the center point of a pixel or voxels of the image and ${\bf{{c}}}_i$ the color coding at the center point.

As in the HK model we define the neighborhood of the particles by

for $i = 1,\dots,n$, and $A_{ij}(t,\epsilon_1,\epsilon_2)$ are entries of the time–varying matrix ${\bf{{A}}}(t,\epsilon_1,\epsilon_2)\in \mathbb{R}^{n\times n}$ defined as

with $\sigma_i$ defined analogously to (4). Here, $\epsilon_1\geq 0$ and $\epsilon_2\geq 0$ are two bounded confidence levels. The two metrics $\left\lVert{\cdot}\right\rVert_{ \mathbb{R}^{d_1}}$ and $\left\lVert{\cdot}\right\rVert_{ \mathbb{R}^{d_2}}$ need to be properly defined according to the specific context of the problem. Then, the mathematical model for any $t\geq 0$ is given by

and initial condition $c_{i,k}(0) = c_{i,k}^0$ and $x_{i,k}(0) = x_{i,k}^0.$ Notice that a sufficient condition to reduce model (8) to the multi–dimensional version of the HK model (1) is $\epsilon_2 > \max_{i,j = 1,\dots,n} \left\lVert{{\bf{{c}}}_i-{\bf{{c}}}_j}\right\rVert_{ \mathbb{R}^{d_2}}$.

2.1.   Properties of the microscopic model

Existence and convergence of solutions to system (8) can be established by using same techniques as in [8] where the original HK in the one–dimensional case was analyzed. In fact, system (7)- (8) belongs to the same class of state-switched systems.

In the case of symmetric interactions we can recover from the previous model similar results on the moment behavior as presented for the HK model in [8,26]. We only record the results here since the proofs are slight variations of existing results (see for example [8]).

Define the moments ${\bf{{m}}}_1(t) \in \mathbb{R}^{d_1}$ and ${\bf{{m}}}_2(t) \in \mathbb{R}^{d_1\times d_1}$ with respect to the time dependent feature as

then the following result holds true.

Lemma 2.1. Let $({\bf{{x}}}_i(t))_{1\leq i \leq n}$ be the solution of the dynamical system (8) with symmetric interactions, i.e. $\sigma_i = n$ in (6). Then, we obtain

Corollary 1. Let $({\bf{{x}}}_i(t))_{1\leq i \leq n}$ be the solution of the dynamical system (8). Assume that $\epsilon_1$ and $\epsilon_2$ are sufficiently large so that interactions are global. Then the first moment is conserved and the following decay estimate holds:

Later, we will show that similar results hold for the continuous model. Some remarks on further properties concerning the particle model (8) are in order.

Remark 1. ● Formation of clusters in the large time behavior is extensively investigated in the review article [35] for general systems of the form (8). However, number of clusters cannot be a–priori predicted starting from a given initial configuration.

● In the non–symmetric case, conservation of the first moment does not hold true. Also, in general, it is not possible to show the decay of the second moment although clustering still appears, see [35].

● Extensions of the model to take into account non static features are obtained by including an interaction term on the right hand side of (7). The determination of such interaction term, however, would be rather problem dependent and in this paper we will not explore further this direction.

3.   Mean–field description

In the case of many particles we derive the formal mean–field equation. Let $\Omega_1 \subseteq \mathbb{R}^{d_1}$, $\Omega_2 \subseteq \mathbb{R}^{d_2}$ be compact domains and $\Omega = \Omega_1 \times \Omega_2$. For $n\geq 1$ we denote by $f_n:\mathbb{R}^+\times \Omega \to \mathbb{R}$ the empirical distribution on $\Omega \subset \mathbb{R}^{d_1\times d_2}$ given by

Let us consider a test function $\varphi({\bf{{x}}},{\bf{{c}}}) \in C^1_0(\Omega)$, i.e. the space of continuous and compactly supported functions on $\Omega$ with continuous derivative. Denote by $\langle \cdot,\cdot \rangle$ the integration of $f_n$ against the test function $\varphi$ on $\Omega$. We have

where we defined

and we used the fact that equation (8) can be re-written as

with $\sigma_i = n \sigma(t,{\bf{{x}}}_i(t),{\bf{{c}}}_i(t))$ and

We can easily compute

and similarly

Collecting these formal computations, after integration by part in ${\bf{{x}}}$, we obtain the weak form of the mean–field equation

If we now define a kernel $\mathcal{A}$ as continuous extension of the adjacency matrix (6)

and ${\bf{{V}}}$ given by

in the limit $n\to\infty$, assuming that the empirical measure $f_n(t,{\bf{{x}}},{\bf{{c}}})$ converge to $f(t,{\bf{{x}}},{\bf{{c}}})$, we formally obtain the strong form of the kinetic equation as

Rigorous analytical results on convergence in the case of $\epsilon_2$ very large and symmetric interactions have been already discussed, for instance in [13,14]. These results guarantee convergence of the distribution $f$ in (12) to a probability limit distribution $f^\infty$, provided the initial distribution $f_0$ at time $t = 0$ has finite second moment and $\mathcal{A}$ is a non–negative, bounded, measurable and symmetric kernel. Observe that these assumptions are verified also in our framework. For a more detailed discussion on analytical results on convergence in mean–field models, we refer e.g. to [12,16]

We briefly discuss the relation to similar kinetic models.

In [9] the analysis of a homogeneous kinetic model for opinion dynamics under bounded confidence is studied. Therein, the model is derived by using a Boltzmann–like approach with instantaneous binary interactions describing compromise. An analogous derivation of the kinetic equation (12) is also possible using a binary interaction model based on an explicit Euler discretization of the underlying particle dynamics (8) with $n = 2$ and performing a grazing collision limit. We omit this computation.

In [9] the authors prove the weak convergence of the solution to a convex combination of Dirac delta functions. In [13] a similar asymptotic distribution is found for the mean–field limit of the classical Hegselmann–Krause model.

A further symmetric clustering model with weighted interactions with respect a fixed number of closest neighbors and corresponding mean–field limit has been introduced in [3]. Moments and long–time behavior could also be studied therein. Finally, models for other applications are also able to cluster information e.g. in traffic flow modeling [42] where the physical acceleration has the role of the bounded confidence.

3.1.   Properties of the mean–field model

As preliminary remark we observe that the marginal distribution $\tilde{f}^c(t,{\bf{{c}}}): = \int_{\Omega_1} f(t,{\bf{{x}}},{\bf{{c}}}) \mathrm{d}{\bf{{x}}}$ is preserved in time by the kinetic equation (12). Instead, it is easy to check that the marginal distribution $\tilde{f}^x(t,{\bf{{x}}}): = \int_{\Omega_2} f(t,{\bf{{x}}},{\bf{{c}}}) \mathrm{d}{\bf{{c}}}$ is not stationary and its behavior in time is still influenced by a ${\bf{{c}}}$ dependent kernel. These considerations are direct consequence of the microscopic model (8) and (7). For this reason, in the following results and if not otherwise stated, we mainly focus on the analysis of moments with respect to the variable ${\bf{{x}}}$.

As in the discrete case we define the $p$–th moment of the kinetic distribution with respect to ${\bf{{x}}}$ as

Here, we used the following multi–index notation ${\bf{{x}}}^{\boldsymbol{{\alpha}}} = x_1^{\alpha_1} \cdots x_{d_1}^{\alpha_{d_1}}$ and $\left|{\boldsymbol{{\alpha}}}\right| = \sum_{i = 1}^{d_1} \alpha_i$. In particular, for each $k,j = 1,\dots,d_1$ we will denote

the first moment and the second moment, respectively. Then, for a symmetric kernel $\mathcal{A}$ we obtain by the integration of the kinetic equation

Lemma 3.1. Let $f(t,{\bf{{x}}},{\bf{{c}}})$ be the solution of the model (12) with symmetric kernel $\mathcal{A}$, i.e., $\sigma(t,{\bf{{x}}},{\bf{{c}}}) = 1$. Then, the following a–priori estimates hold true:

for each $k,j = 1,\dots,d_1$.

Proof. For each fixed $k = 1,\dots,d_1$ we have

The last term vanishes due to anti–symmetry of the integrand by the interchange of variables ${\bf{{x}}} \leftrightarrow {\bf{{y}}}$. The second statement follows by observing that for each $k = 1,\dots,d_1$ we have

Hence, we obtain

Instead, in the case of global interactions and for a general kernel $\mathcal{A}$ and we have

Corollary 2. Assume that $\epsilon_1$ and $\epsilon_2$ are sufficiently large so that $\mathcal{A}(t,{\bf{{x}}},{\bf{{c}}},{\bf{{y}}},{\bf{{z}}};\epsilon_1, \epsilon_2) = 1$, for all ${\bf{{x}}},{\bf{{y}}}\in\Omega_1$, ${\bf{{c}}},{\bf{{z}}}\in\Omega_2$ and $t\geq 0$. Then, for all $k,j = 1,\dots,d_1$, the following relations hold true:

Proof. Multiplying (12) by ${\bf{{x}}}^{\boldsymbol{{\alpha}}}$ with $\left|{\boldsymbol{{\alpha}}}\right| = p\in\mathbb{N}$, and integrating over $\Omega$, we compute

where $\boldsymbol{{\alpha}}^{(\ell)} = \boldsymbol{{\alpha}} - \bf{{e}}_\ell$ with $\bf{{e}}_\ell$ being the $\ell$–th element of the standard basis vector of $\mathbb{R}^{d_1}$ and $\left|{\boldsymbol{{\alpha}}^{(\ell)}}\right| = p-1$. Then, for $p = 1$ we still have have conservation of the first moment since we obtain

For $p = 2$ we have

and therefore

Observe that, compared to the local interaction case analyzed in Theorem 3.1 where in general the energy decay property of all second order moments is not guaranteed, in global interactions also the mixed second order moments decay in time. In other words, while in the local interaction model only variances go to zero, in the global interaction model both variances and covariances go to zero in the large time behavior.

Lemma 3.1 and Corollary 2 suggest that for $t\to \infty$ any initial distribution will tend to a stationary distribution that is concentrated on a finite number of points in $\Omega_1$ and represents therefore clusters. Hence, in the following we investigate the asymptotic distribution of the model (12).

Lemma 3.2. Let $\epsilon_1$ and $\epsilon_2$ be arbitrary positive bounded confidence levels. Then, the distribution

with $f_k>0$ and $\sum_{k = 1}^{n_1} \sum_{\ell = 1}^{n_2(k)} f_k = 1$ is a weak stationary solution of the model (12) if and only if either $\left\lVert{{\bf{{x}}}_i-{\bf{{x}}}_k}\right\rVert_{ \mathbb{R}^{d_1}}>\epsilon_1$ for all $i\neq k$ or $\left\lVert{{\bf{{c}}}_{i_j}-{\bf{{c}}}_{k_\ell}}\right\rVert_{ \mathbb{R}^{d_2}}>\epsilon_2$ for all $i\neq k$, for all $j,\ell$, or both hold true.

Proof. Let us assume that at least one of $\left\lVert{{\bf{{x}}}_i-{\bf{{x}}}_k}\right\rVert_{ \mathbb{R}^{d_1}}>\epsilon_1$ for all $i\neq k$ and $\left\lVert{{\bf{{c}}}_{i_j}-{\bf{{c}}}_{k_\ell}}\right\rVert_{ \mathbb{R}^{d_2}}>\epsilon_2$ for all $i\neq k$, for all $j,\ell$, is verified and prove that $f^\infty$ is the asymptotic distribution of the equation (12), i.e. given a test function $\varphi\in C_0^1(\Omega)$

We obtain

Conversely, assume that $f^\infty$ as in (13) is the steady state of (12). Assume by contradiction that there exist $k,i\in\{1\dots,n_1\}$ such that $\left\lVert{{\bf{{x}}}_k-{\bf{{x}}}_i}\right\rVert_{ \mathbb{R}^{d_1}}\leq\epsilon_1$ and $\ell\in\{1\dots,n_2(k)\}$, $j\in\{1\dots,n_2(i)\}$ such that $\left\lVert{{\bf{{c}}}_{i_j}-{\bf{{c}}}_{k_\ell}}\right\rVert_{ \mathbb{R}^{d_2}}\leq\epsilon_2$. Then, using similar computations as before we obtain

which contradicts the hypothesis.

Lemma 3.3. Assume that $\epsilon_1$ and $\epsilon_2$ are sufficiently large so that $\mathcal{A}(t,{\bf{{x}}},{\bf{{c}}},{\bf{{y}}},{\bf{{z}}};\epsilon_1,\epsilon_2) = 1$, for all ${\bf{{x}}},{\bf{{y}}}\in\Omega_1$, ${\bf{{c}}},{\bf{{z}}}\in\Omega_2$ and $t\geq 0$. Then, the distribution

with $\tilde{f}^c(0,{\bf{{c}}}) = \int_{\Omega_1} f_0({\bf{{x}}},{\bf{{c}}}) \mathrm{d}{\bf{{x}}}$ is a weak stationary solution of the model (12) if and only if $\bar{x}_k = u_k(0)$.

Proof. Substituting the expression of $f^\infty$ in the weak form of the kinetic equation, integrating by parts and using the fact that for $\epsilon_1$ and $\epsilon_2$ large enough so that

we obtain

which is zero if $\bf{{\bar{x}}} = \int_{\Omega} {\bf{{y}}} f_0({\bf{{y}}},{\bf{{z}}}) \mathrm{d}{\bf{{z}}} \mathrm{d}{\bf{{y}}}$.

Remark 2. As direct consequence of Theorem 3.2 and Theorem 3.3 we have that taking $\epsilon_2$ sufficiently large or an initial distribution $f_0({\bf{{x}}},{\bf{{c}}})$ being atomic with respect the second variable, the same quantized steady–state of the kinetic formulation of the microscopic model (1) are preserved, cf. [13].

Remark 3. In Theorem 3.2 the number of clusters $n_1$, as well the values $f_j$ and the positions ${\bf{{x}}}_j$ of the clusters, are functions of the initial distribution $f_0({\bf{{x}}},{\bf{{c}}})$ and of the bounded confidence levels $\epsilon_1$ and $\epsilon_2$. As pointed–out also in [14], in general it is not possible to predict the number of Dirac deltas in the asymptotic configuration from a given initial distribution. However, for $\epsilon_2$ sufficiently large and assuming $\Omega = [0,1]^d$, the number of clusters is $1\leq \tilde{n} \leq \lfloor \frac{1}{\epsilon_1^d} \rfloor$. This consideration is also observed at the Boltzmann level, see [9].

4.   Numerical experiments and applications

The theoretical results on the asymptotic behavior of the mean–field equation (12) introduced in the previous section are here also numerically investigated. Moreover, we show the efficiency of the model as technique to solve realistic data–clustering problems and to this end we focus on applications to the field of shape detection and image segmentation.

In order to efficiently solve the kinetic model (12) we employ the Mean Field Interaction Algorithm introduced in [2] which is based on random subset evaluation of the kernel term $\mathcal{A}_{\epsilon_1,\epsilon_2}(t,{\bf{{x}}},{\bf{{c}}},{\bf{{y}}},{\bf{{z}}})$ given in (10). The algorithm used in the numerical experiments is summarized by the steps in Algorithm 1 for a time interval $[0,T]$ discretized in $k_{tot}$ subintervals of size $\Delta t$. The computational cost of the algorithm is $O = (MN)$, where $M$ is the size of the subset of interacting particles, and for $M = N$ we obtain the explicit Euler scheme for the original $N$ particle system (8) whose cost is $O(N^2)$. For further details on the Mean Field Interaction Algorithm we refer to [2].

4.1.   Numerical steady–states and moment evolution

We use the Mean Field Interaction Algorithm to numerically investigate the properties of the kinetic model proved in the previous section. We analyze two typical situations which lead to the two applications we show later in this section.

4.1.1.   Constant static feature.

First, we consider an initial distribution being the uniform distribution with respect to ${\bf{{x}}}$ and a constant distribution along the static feature variable ${\bf{{c}}}$, i.e. $f_0({\bf{{x}}},{\bf{{c}}}) = \chi_{[0,1]^{d_1}}(\bf{{\vec{x}}})$. This choice is obviously also consistent with considering the bounded confidence level $\epsilon_2$ very large. We show that, as discussed in Remark 2, equation (12) provides the same steady–states of the Hegselmann–Krause model.

One–dimension. The numerical steady–states of the mean–field equation (12) are first computed for the one–dimensional case and compared to the steady–states of the microscopic model (8). Observe that, under the assumption of an initial constant distribution along ${\bf{{c}}}$, the asymptotic behavior of (8) is the classical one prescribed by the one–dimensional Hegselmann–Krause model (1). The evolution in time of the first and second moment is also provided.

In Figure 1 we show the steady–state provided by the mean–field kinetic model (12) for the one–dimensional initial uniform distribution on $[0,1]$. The final time is $T = 20$ and the time step is $\Delta t = 0.5$. We consider $N = 5\times 10^5$ particles in the Monte Carlo method so that we reduce error due to the sampling. The number of interacting particles is taken as $M = 10$. We show the results for two values of the confidence bound, $\epsilon_1 = 0.5$ in the left panel and $\epsilon_1 = 0.15$ in the right panel. The evolution of the distributions up to final time is showed by normalizing with respect to the maximum value at each fixed time. For $\epsilon_1 = 0.5$ we observe the formation of a consensus state in large time behavior. In fact, the final distribution is a Dirac delta centered in the initial value of the first moment. For $\epsilon_1 = 0.15$, instead, three clusters arise at equilibrium, similarly the classical Hegselmann–Krause model for the same confidence bounds and equally distributed initial data.

As observed in Figure 2, for both values of the bounded confidence level, we have conservation of the first moment and energy dissipation in time. Although we are using a non–symmetric model, the first moment is preserved during the time evolution, with deviation of order $10^{-4}$ from the initial value, since the initial distribution is symmetric. The fluctuations in the first moment are due to the stochastic method.

In Figure 3 we show some results of a simple analysis of Algorithm 1 with respect to the values of $M$, i.e. the size of the random subset of interacting particles. Left panel shows the trend to the steady–state of the kinetic model computed with Algorithm 1 and same parameters as in the right plot of Figure 1. However, here we consider $N = 2\times10^4$ and $M = 2$. The same equilibrium state is reached but with a faster transient. The dependence of the velocity to the formation of cluster on the size $M$ of the random subset is showed in the right panel of Figure 3. The velocity to equilibrium decreases as $M$ increases. In particular, the $M = 2$ case, where each particle is enforced to align to the velocity of the other particle instead of their "mean", exhibits the fastest convergence. This analysis suggests that, although model (12) has a convolution structure and in principle can be efficiently solved with a Fast Fourier Transform at a $O(N\log N)$ cost, Algorithm 1 may be preferable since it has a comparable (or even lower) computational cost and it can be applied to more general kernels where the convolution structure is lost.

Two–dimensions. We also show one example of the numerical steady–states of the mean–field equation (12) in the two–dimensional case. Again we employ Algorithm 1. Figure 4 shows the steady–state for the particle density and the kinetic density at time $t = 4$ and final time $T = 50$. We use $N = 10000$ particles and the time step is $\Delta t = 0.5$. Using the Mean Field Interaction Algorithm 1, the number of interacting particles is taken as $M = 10$. We show the results for one value of the confidence bound, $\epsilon_1 = 0.15$ which leads to the formation of $8$ clusters at equilibrium. The evolution in time of the two–dimensional moments is also provided in Figure 5. We still observe conservation of the first moment, due to the initial symmetric distribution, and decay in time for the two non–mixed second moments.

4.1.2.   Non–constant static feature.

Next, we consider a non–constant initial distribution along the variable ${\bf{{c}}}$. We consider a one–dimensional setting both along ${\bf{{x}}}$ and ${\bf{{c}}}$ since the goal of this experiment is to provide evidence of the influence of the static feature on the clustering process. The initial kinetic distribution it given as a tensor product between a uniform and Gaussian distribution so that

with mean $\mu = 0.5$, and variance $\sigma^2 = 0.3$. See the top left plot in Figure 6. All the simulations are performed by using $N = 5000$ samples up to equilibrium.

In this experiment, we expect that the behavior at equilibrium is influenced also by the interactions with respect to the static feature variable due to the corresponding characteristic function in the kernel (10). In Figure 6 we show the results provided by using $\epsilon_1 = 0.15$ and $\epsilon_2 = 0.025$. The top left plot shows particles and kinetic density at equilibrium. Observe that, compared to the case with constant distribution along the static feature, the low value of $\epsilon_2$ allows for the formation of more clusters, precisely $8$. Clusters arise because either their distance in $x$ direction is larger than the confidence level $\epsilon_1$ or the minimum distance in $c$ direction is larger than $\epsilon_2$. This consideration is analyzed in the bottom left plot of Figure 6 where the blue line with circle markers shows the distance in $x$ between two adjacent clusters and the red line with triangle markers the minimum distance between the static features of particles being in two adjacent clusters. Observe that when the distance in $x$ is lower than the value $\epsilon_1$, the corresponding distance between the static features is larger than $\epsilon_2$, and therefore cluster is no longer possible. Moreover, we notice that time to reach equilibrium is larger due to two levels of clustering. Finally, the bottom right plot of Figure 6 shows the two marginal distributions.

In Figure 7 we show the behavior at equilibrium by considering two different pairs of confidence levels. In the left plot, $\epsilon_1 = 0.15$ and $\epsilon_2 = 0.1$. Here $\epsilon_2$ is large enough to not influence the clustering process and in fact we recover $3$ clusters, exactly the same provided in Figure 1 when the initial kinetic distribution was constant in $c$. In the right plot, $\epsilon_1 = 1$ and $\epsilon_2 = 0.025$. In this case, $\epsilon_1$ is very large and we show the formation of clusters at equilibrium when only clustering due to the static feature is allowed.

4.2.   Clustering and shape detection

The property of the equation (12), analytically studied in Section 3.1 and numerically investigated in Section 4.1, of having quantized steady–states makes the model suitable for solving data clustering problems. In particular, in this section we use the kinetic model as technique for shape detection.

Shape detection can be considered as a branch of pattern recognition problems that focuses on the discovering of previously unknown patterns in a big set of data. Machine learning and many clustering algorithms have already been applied to this class of problems, such as the $k$–means algorithm which however suffers from the necessity of defining a–priori the number of clusters. Instead, as already pointed out, opinion dynamics based models automatically find the number of clusters as functions of the confidence value.

More specifically, in this section we apply the models introduced in Section 2 to the example of a character recognition. We consider a letter "A" which is composed by $n$ two-dimensional points $\mathcal{L} = \{{\bf{{x}}}_i\}_{i = 1}^n$ defining three segments in $\Omega = [0,1]^2$. To each of these points we apply an additive noise distributed according to a uniform or a Gaussian distribution. We obtain a new set of $n$ points

which define the noisy pattern. Here, the value $\alpha$ represents the percentage of noise and it is also chosen in such a way that $\tilde{{\bf{{x}}}}_i\in\Omega$, $i = 1,\dots,n$. The points $\{\tilde{{\bf{{x}}}}_i\}_{i = 1}^n$ are used as initial condition of the model. They are seen as samples from a distribution defining the noisy pattern and the Mean Field Interaction Algorithm 1 is employed in order to solve the kinetic model.

The goal of this example is to cluster the noisy data into a set of points which give information on the shape of the exact pattern to be detected. This can be considered also as a dimensionality reduction problem which could help other algorithms to recognize the unknown pattern efficiently by using information on the position of the clusters. Here the efficiency of the results is studied by means of the following measure:

where $\{ {\bf{{C}}}_k \}_{k = 1}^{\tilde{n}}$ is the position of clusters in $\Omega$, $\tilde{n}$ the number of clusters at equilibrium. The quantity $\mathcal{E}$ measures for each cluster the minimum $2$–norm distance to the exact pattern and then computes the normalized $1$–norm of these quantities. The definition of the measure $\mathcal{E}$ is inspired by the idea of computing an average of the minimum distances between the clusters and the shape $\mathcal{L}$, when formation of multiple clusters distributed close the shape $\mathcal{L}$ arises, as in the cases showed here. Certainly, $\mathcal{E}$ cannot provide a general way to measure the quality of the results, since e.g. in case of one steady-state cluster positioned exactly on a point of $\mathcal{L}$, we would have $\mathcal{E} = 0$ providing absurdly the optimal choice.

In the following, all the results are given for an initial set of $n = 5000$ particles and a large enough final time to reach a steady–state $T = 50$. Moreover, we use the Euclidean norm and all particles are initialized with a constant static feature. For an extension of this application to clustering also with respect to the static feature see Remark 4 below.

In Table 1 and Table 2 we show a sensitivity study of the error $\mathcal{E}$ as function of the percentage of noise $\alpha$ and of the bounded confidence value $\epsilon_1$, for initial noisy data uniformly and normally distributed, respectively. The number of clusters at equilibrium is also provided. Results with minimal error value are highlighted with bold font. In particular, we point out as the increasing percentage of noise highly affects the results.

In Figure 8 and Figure 9 some of this cases are shown. Precisely, for the uniformly distributed initial noisy data, Figure 8, we consider the largest value of the percentage of noise $\alpha = 10\%$ and three values of the bounded confidence level $\epsilon_1$: the smallest $\epsilon_1 = 0.06$ in Table 1, the optimal $\epsilon_1 = 0.08$ (giving the smallest value of the error $\mathcal{E}$) and the largest one $\epsilon_1 = 0.1$. We observe that the result is highly influenced by the choice of the confidence bound which determines the number of final clusters. Smaller and larger values of $\epsilon_1$ result in larger values of the error measure $\mathcal{E}$. In fact, in one case the model provides too many clusters which are very spread out, with several outliers. In the other case too few clusters so that all the information on the exact shape is lost. For the optimal value of $\epsilon_1 =$ the model provides a good agreement with the exact pattern.

Similar considerations hold for Figure 9 where we show the noisy initial data obtained with $\alpha = 10\%$ percentage of noise and three values of the bounded confidence level $\epsilon_1 = 0.05$, $\epsilon_1 = 0.0675$ (giving the smallest value of the error $\mathcal{E}$) and $\epsilon_1 = 0.08$.

We point out that the goal of this application is to provide a preliminary step for shape feature extraction. This technique should be coupled with a further algorithm which is able to select the correct letter from a given alphabet. The dimensionality reduction provided by the data clustering approach allows to efficiently employ the subsequent recognition analysis.

Remark 4 (Shape detection with non–constant static feature). In the previous examples the static features is taken constant for all the data. In the non–constant case, one could think of this feature as an additional initial given input which measures the quality of the information of each single point and that can be modeled by the distance to its exact value in $\mathcal{L}$.

4.3.   Color image segmentation

We now turn to the second application of the data clustering model 12 which is the segmentation of gray scale images. Image segmentation is widely used in medical and astronomical image processing, face recognition, etc. In fact, this technique allows to partition an image in sets of significant regions of pixels sharing same characteristics such as closeness and similar intensity of color. As in shape detection, the aim of such approach is to modify the description of an image into a structure which easier to be analyzed for subsequent computer vision algorithms.

Several mathematical techniques have been applied to image segmentation. For instance we mention level set methods [39,43] methods based on the Kuramoto model [29,37] and supervised convolution neural networks [17,46]. A more complicated procedure to segmentation of images is the clustering approach and we refer e.g. to the $k$–means method [44,47], $c$–means method [33] and hierarchical clustering method [45].

Here, we study the efficiency of the kinetic model (12) to solve image segmentation problems and therefore we still relies on the clustering technique. Thanks to the model introduced in this paper, segmentation can be performed by using two levels of clustering in order to determine regions of pixels with similar characteristics. More precisely, we cluster based on the Euclidean distance and the distance of the gray intensity between pixels. Taking into account spatial coordinates is needed in order to avoid the problem of selecting homogeneous regions of pixels which are however distinct in the original image. The clustering with respect the intensity of the gray colors is performed by using the static feature variable ${\bf{{c}}}$.

The application of (12) works as described in the following. Number of pixels is the number of particles which are assumed to be equally spaced samples. The intensity of the gray color of each pixel is computed at initial time and defines the one–dimensional static feature $c$ of each particle. The Mean Field Interaction Algorithm 1 is then applied to find clusters representing local regions of the original image with homogeneous intensity of color. At equilibrium, we compute the mean of the color intensity of pixels belonging to the same cluster and then they are mapped back to the original positions in the initial image to get the segmentation. Another possible approach, which is more suitable for images with very sharp color intensities, such as in astronomical images, is to apply thresholding technique. This method is based on replacing the color intensity of each cluster by black if the mean is below a certain threshold or by white otherwise. This technique results therefore in binary images.

In Figure 10 we first apply the segmentation to a benchmark test in order to show how the model works on a simple test. We initialize an gray scale image with $4096$ pixels with four regions of different intensity colors. The top and the bottom left corners have intensity $1$ and $0.75$, respectively. The top and the bottom right corners have intensity $0$ and $0.25$, respectively. See the left panel of Figure 10. The positions of the pixels are rescaled on $\Omega = [0,1]^2$. The clustering algorithm is applied with the two bounded confidence values $\epsilon_1 = 0.5$ and $\epsilon_2 = 0.3$. At equilibrium the two red clusters showed in the middle panel of Figure 10 are computed by the kinetic model. The mean of the intensity colors of the pixels in each cluster is computed and mapped back to the initial position obtaining the segmentation in the right panel of Figure 10. The two values of the gray intensity are $0.125$ and $0.875$.

We now apply the segmentation process to real images taken from different data–sets. In Figure 11 we consider an image of $174\times73$ pixels showing a night sky, trees and the moon. On this example, the method is able to provide good results with a wide range of confidence values since the variation of the gray scale is sharp between the objects which are also characterized by regions having homogeneous colors. In Figure 11 we report the result obtained with $\epsilon_1 = 0.2$ and $\epsilon_2 = 0.1$ which results in the formation of $5$ clusters.

The second image in Figure 12 is also characterized by two backgrounds: the sky which is very homogeneous and the field which is less homogeneous. The goal of the segmentation process would be to identify the two cows, separating them from the background. We report the results obtained with $\epsilon_1 = 0.2$ and $\epsilon_2 = 0.05$, resulting in $6$ clusters. Observe that $\epsilon_2$ is taken small in order to avoid the formation of many clusters due to the low homogeneity of the field. The two cows are well identified by the method.

In the case of Figure 13 the confidence level $\epsilon_2$ is taken larger in order to force clustering of regions with snow but different shadows. The clustering method identifies $8$ clusters and the skier is well separated by the background.

In Figure 14 we employ the kinetic model for clustering to segmentation of a real image taken by [32]. The true image is $67\times67$ pixels and several values of the bounded confidence levels $\epsilon_1$ and $\epsilon_2$ are applied in order to cluster with respect positions and intensity color of pixels. The best result is obtained with $\epsilon_1 = 0.3$ and $\epsilon_2 = 0.1$ which results in the formation of $8$ clusters. Lower values of $\epsilon_1$ allows to obtain more clusters and thus more gray intensity colors at equilibrium. Increasing $\epsilon_1$ requires to not take $\epsilon_2$ also large, otherwise many information are lost at equilibrium.

Finally, in Figure 15 we present an additional experiment aimed to provide a comparison between our segmentation process and the one published in [29] based on the Kuramoto model. The image is selected from the data–set [32]. Right panel of Figure 15 has to be compared with Figure (7c) in [29] and, in particular, we observe that our result is able to detect better some structures of the image that are lost by using the Kuramoto model.

5.   Conclusions

In this paper we have proposed a new method for clustering of large data which is based on a suitable generalization of the Hegselmann-Krause opinion dynamics model. The extension accounts for clustering with respect to two different sets of features of each single datum: one set describes time dependent characteristics, the other one represents static characteristics.

The mean–field limit of the particle model has been formally derived and the resulting kinetic equation allows for the study of mathematical properties. Time evolution of moments and asymptotic behavior of kinetic equation have been analytically characterized. Moreover, the derivation of the mean–field model allows for reduction of computational complexity thanks to a suitable random subset Monte Carlo algorithm.

Several applications to digital imaging have been proposed. In particular, we have focused on application to shape detection and image segmentation, showing the efficiency of the model to provide satisfying results. Finally, we emphasize that the present model has to be intended as a starting point towards more realistic applications, for example based on the use of non static features combined with suitable machine learning approaches.

Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy – EXC-2023 Internet of Production – 390621612.

Lorenzo Pareschi and Giuseppe Visconti acknowledge the support of the "National Group for Scientific Computation (GNCS-INDAM)", project "Numerical approximation of hyperbolic problems and applications".