Robust Clustering on High-Dimensional Data with Stochastic Quantization
Anton Kozyriev Vladimir Norkin

    

    
    Anton Kozyriev and Vladimir Norkin

    National Technical University of Ukraine ”Igor Sikorsky Kyiv Polytechnic Institute”, Kyiv, 03056, Ukraine a.kozyriev@kpi.ua V.M.Glushkov Institute of Cybernetics, National Academy of Sciences of Ukraine, Kyiv, 03178, Ukraine v.norkin@kpi.ua

    Robust Clustering on High-Dimensional Data with Stochastic Quantization
    Anton Kozyriev1 Vladimir Norkin1,2
    September 9, 2024
=======================================================================



    § ABSTRACT 

This paper addresses the limitations of traditional vector quantization (clustering) algorithms, particularly K-Means and its variant K-Means++, and explores the Stochastic Quantization (SQ) algorithm as a scalable alternative for high-dimensional unsupervised and semi-supervised learning problems. Some traditional clustering algorithms suffer from inefficient memory utilization during computation, necessitating the loading of all data samples into memory, which becomes impractical for large-scale datasets. While variants such as Mini-Batch K-Means partially mitigate this issue by reducing memory usage, they lack robust theoretical convergence guarantees due to the non-convex nature of clustering problems. In contrast, the Stochastic Quantization algorithm provides strong theoretical convergence guarantees, making it a robust alternative for clustering tasks. We demonstrate the computational efficiency and rapid convergence of the algorithm on an image classification problem with partially labeled data, comparing model accuracy across various ratios of labeled to unlabeled data. To address the challenge of high dimensionality, we trained Triplet Network to encode images into low-dimensional representations in a latent space, which serve as a basis for comparing the efficiency of both the Stochastic Quantization algorithm and traditional quantization algorithms. Furthermore, we enhance the algorithm's convergence speed by introducing modifications with an adaptive learning rate.




    

§ INTRODUCTION


Quantization and clustering are fundamental encoding techniques that provide compact representations of original data <cit.>. Clustering algorithms have emerged as prominent tools for unsupervised learning, with recent applications spanning diverse domains such as location-allocation problems <cit.>, document classification <cit.>, and data compression <cit.>.

In this context, we consider a random variable ξ with values in Euclidean space ℝ^n and distribution P(dξ), representing the original distribution. The encoded discrete distribution is parameterized by a set of atoms {y_1, …, y_K} with corresponding probabilities {q_1, …, q_K}. The optimal quantization problem aims to find the encoded distribution that minimizes the distance to the original distribution. This mathematical structure is analogous to the optimal clustering problem, where the objective is to determine the positions of K cluster centers {y_1, …, y_K} such that the sum of distances from each element ξ to the nearest cluster center is minimized.

The K-Means algorithm, proposed by Lloyd <cit.>, has been widely used for solving quantization and clustering problems, with numerous extensions <cit.>. Bottou and Bengio <cit.> interpreted the K-Means algorithm as an analogue of Newton's method and proposed several stochastic gradient descent algorithms for solving optimal quantization and clustering problems. However, traditional clustering algorithms are limited by the requirement to load all training data into memory, rendering them non-scalable for large datasets. To address this limitation, Sculley <cit.> introduced the Mini-Batch K-Means algorithm, which utilizes only a small subset of possible ξ values at each iteration.

The Stochastic Quantization algorithm reframes the clustering problem as a stochastic transportation problem <cit.> by minimizing the distance between elements of the original distribution {ξ} and atoms of the encoded discrete distribution {y_k}. This approach employs Stochastic Gradient Descent (SGD) <cit.> to search for an optimal minimum, leveraging its computational efficiency in large-scale machine learning problems <cit.>. The use of stochastic approximation allows the algorithm to update parameters with only one element ξ per iteration, ensuring memory efficiency without compromising convergence to the minimum <cit.>.

This paper explores advanced modifications of the Stochastic Quantization algorithm, incorporating accelerated variants of SGD <cit.> and adaptive learning rate techniques <cit.> to enhance convergence speed. Norkin et al. <cit.> provide a comprehensive comparison of various SGD variants, highlighting their respective advantages and limitations, while also offering convergence speed estimations.

Given that the optimal quantization problem is non-smooth and non-convex, specialized methods are required <cit.>. To validate its convergence, we apply the general theory of non-smooth, non-convex stochastic optimization <cit.>. While traditional clustering algorithms lack theoretical foundations for convergence guarantees and rely primarily on seeding techniques <cit.>, stochastic optimization theory provides specific conditions for the local convergence of the Stochastic Quantization algorithm, which we supplement in this research.

This paper introduces a novel approach to address semi-supervised learning challenges on high-dimensional data by integrating the Stochastic Quantization algorithm with a deep learning model based on the Triplet Network architecture <cit.>. The proposed method encodes images into low-dimensional representations in the ℝ^3 latent space, generating meaningful encoded features for the Stochastic Quantization algorithm. By employing the Triplet Network, this approach overcomes the limitations of quantization and clustering algorithms in high-dimensional spaces, such as visualization difficulties and decreased precision as the number of dimensions increases <cit.>.

To illustrate the efficiency and scalability of the Stochastic Quantization algorithm, we conducted experiments on a semi-supervised image classification problem using partially labeled data from the MNIST dataset <cit.>. The Triplet Network is initially trained on the labeled portion of the dataset as a supervised learning model. Subsequently, the trained network is utilized to project the remaining unlabeled data onto the latent space, which serves as input for training the Stochastic Quantization algorithm. The performance of the proposed solution is evaluated using the F1-score metric <cit.> for multi-label classification across various ratios of labeled to unlabeled data.

    

§ STOCHASTIC QUANTIZATION


Unlike traditional clustering methods that minimize the distance between each element of Ξ = {ξ_i, i = 1, …, I} and the nearest center Y = {y_k, k = 1, …, K}, Stochastic Quantization conceptualizes the feature set Ξ and cluster centers Y as discrete probability distributions. The Wasserstein (or Kantorovich–Rubinstein) distance is employed to minimize distortion between these distributions when representing a continuous distribution by a discrete one <cit.>. Subsequent research <cit.> has explored the application of quantization algorithms to solve optimal allocation problems for service centers, where each atom of the discrete distribution represents the location of facilities and customers, respectively.

     <cit.>. Optimal quantization minimizes the weighted sum of distances between elements of the feature set {ξ_i}⊂ℝ^n and centers {y_k}⊂ℝ^n:
    
    min_y = { y_1, …, y_K }∈ Y^K ⊂ℝ^nKmin_q = { q_1, …, q_K }∈ℝ^K_+min_x = { x_ij≥ 0 }∑_i=1^I ∑_k=1^K d(ξ_i, y_k)^r x_ik

    subject to constraints:
    
    ∑_k=1^K x_ik = p_i,   ∑_k=1^K q_k = 1,    i = 1, …, I

    where p_i > 0, ∑_i=1^I p_i = 1 are normalized supply volumes, x_ik are transportation volumes, d(ξ_i, y_k)_p = ξ_i - y_k _p = (∑_j=1^n | ξ_ij - y_kj |^p)^1/p is the l_p norm defining the distance between elements in the objective function (<ref>), Y ⊂ℝ^n is a common constraint set for variables {y_k, k = 1, …, K}, and n, I, K ∈ℕ.

In this research, we employ the Euclidean norm (p = 2) as the distance metric, defined as d(ξ_i, y_k)_2 = √(∑_j=1^n | ξ_ij - y_kj |^2). The choice of distance metric may vary depending on the problem domain. For instance, the cosine similarity function d(ξ_i, y_j)_cos = cos(ξ_i, y_j) = ξ_i · y_j/ξ_i · y_j is utilized in text similarity tasks <cit.>, while Kolmogorov and Levy metrics are employed for probability and risk theory problems <cit.>.

It is evident that in the optimal plan, all mass at point ξ_i is transported to the nearest point y_k. Consequently, problem (<ref>)-(<ref>) can be reduced to the following non-convex, non-smooth global stochastic optimization problem, with the objective function defined as:

    min_y = { y_1, …, y_K }∈ Y^K ⊂ℝ^nK F(y_1, …, y_k)

where

    F(y) = F(y_1, …, y_k) = ∑_i=1^I p_i min_1 ≤ k ≤ K d(ξ_i, y_k)^r = 𝔼_i ∼ pmin_1 ≤ k ≤ K d(ξ_i, y_k)^r

Here, 𝔼_i ∼ p denotes the expected value over the random index i that takes values {1, …, I} with probabilities {p_1, …, p_I}, respectively.

     In the global optimum y^* = (y_1^*, …, y_K^*) of (<ref>), all {y_1^*, …, y_K^*} belong to the convex hull of elements {ξ_1, …, ξ_I} in the feature set.


    Assume, by contradiction, that there exists some y_k^*^* ∉conv{ξ_1, …, ξ_I}. Consider the projection y̅_k^*^* of y_k^*^* onto conv{ξ_1, …, ξ_I} and points y_k^*^*(t) = (1 - t)y_k^*^* + ty̅_k^*^*, t ∈ [0, 1]. We observe that ∀ξ_i, t ∈ (0, 1]: y_k^*^*(t) - ξ_i < y_k^*^* - ξ_i. If y_k^*^* - ξ_i^* = min_1 ≤ k ≤ Ky_k^* - ξ_i^* for some i^*, then
    
    min{ y_k^*^*(t) - ξ_i^*, min_k ∉ k^* y_k^* - ξ_i^*} < min_k  y_k^* - ξ_i^*

    Thus, y^* = (y_1^*, …, y_K^*) is not a local minimum of the objective function (<ref>). Now, consider the case where y_k^*^* - ξ_i  > min_k  y_k^* - ξ_i for all i. By assumption, min_k  y_k^* - ξ_i' for some i'. The vector y' = (y_1^*, …, y_k^* - 1^*, ξ_i', y_k^* + 1^*, …, y_K^*) satisfies F(y') < F(y^*), contradicting the assumption that y^* is a minimum. This completes the proof. 


For a continuous probability distribution P(dξ), we can interpret the objective function (<ref>) as a mathematical expectation in a stochastic optimization problem <cit.>:

    min_y = { y_1, …, y_K }∈ Y^K ⊂ℝ^nK[F(y_1, …, y_k) = 𝔼 f(y, ξ) = ∫_ξ∈Ξ f(y, ξ) P(d ξ)]

with 

    f(y, ξ) =  min_1 ≤ k ≤ K d(ξ, y_k)^r,

where the random variable ξ may have a multimodal continuous distribution. The empirical approximation of F(y) in (<ref>) is:

    F_N(y) = 1/N∑_i=1^N min_1 ≤ k ≤ K d(ξ_i, y_k)^r

where {ξ_i, i = 1, …, N} are independent, identically distributed initial samples of the random variable ξ. If K = 1, Y is convex, and r≥ 1, then problem (<ref>) is unimodal and reduces to a convex stochastic optimization problem:

    min_y ∈ Y [ F(y) =  𝔼_ĩ∼ p d(ξ_ĩ, y)^r ]

However, for K ≥ 2, the function f(ξ, y) = min_1 ≤ k ≤ K d(ξ, y_k)^r, y = (y_1, …, y_K) is non-smooth and non-convex. In terms of <cit.>, f(ξ, y) is a random generalized differentiable function, its generalized gradient set can be calculated by the chain rule:

    ∂ f(ξ, y)    =   conv{ (0, …, 0, g_k^*, 0, …, 0),    k^* ∈ S(ξ, y),    0 ∈ℝ^n }
    
            S(ξ, y)    =   { k^*: ξ - y_k^* = min_1 ≤ k ≤ Kξ - y_k }
    
            g_k^*   =    r ξ - y_k^*^r - 2 (y_k^* - ξ)


The expected value function (<ref>) is also generalized differentiable, and the set 𝔼_ξ∂ f(ξ, y) is a generalized gradient set of the function F <cit.>. Vectors g(ξ) = (0, …, 0, g_k, 0, …, 0), k ∈ S(ξ, y), 0 ∈ℝ^n, are stochastic generalized gradients of the function F(y_1, …, y_K).


These gradients can be utilized to find the optimal element y_k^* in a feature set Ξ using Stochastic Gradient Descent (SGD) <cit.>:

    y_k+1 = π_Y (y_k - ρ_k g_k^*),   π_Y (x) = _y ∈ Y x - y,    y^0 ∈ Y,    k ∈ℕ,

where ρ_k > 0 is a learning rate parameter, and π_Y is the projection operator onto the set Y. The iterative process (<ref>)-(<ref>) for finding the optimal element is summarized in Algorithm <ref>. While SGD is an efficient local optimization algorithm, the ultimate task is to find global minima of (<ref>). The research in <cit.> proposes a stochastic branch and bound method applicable to the optimization algorithm (<ref>). The idea is to sequentially partition the initial problem into regions (with constraint set Y_1 ×…× Y_K) and use upper and lower bounds to refine partitions with the so-called interchanges relaxation to obtain lower bounds:

    min_{ y_k ∈ Y_k } F(y_1, …, y_K)
               ≥   ∑_i=1^I p_i min_y ∈ Y_1×…× Y_Kmin_1 ≤ k ≤ K d(ξ_i, y_k)^r 
       ≥   ∑_i=1^I p_i min_1 ≤ k ≤ K d(ξ_i, π_k(ξ_i)^r.


The local convergence conditions of the stochastic generalized gradient method for solving problem (<ref>) are determined in Theorem <ref>, with the proof provided in <cit.>.

     <cit.>. Consider the iterative sequence { y^(t) = (y_1^(t), …, y_K^(t)) }:
    
    y_k^(t)   := π_Y (y_k^(t) - ρ_t g_k^(t))     k^(t)   = S(ξ̃^(t), y^(t))     t    = 0, 1, 2, …
    
                g_k^(t)   = r ξ̃^(t) - y_k^(t)^r - 2 (y_k^(t) - ξ̃^(t))     k    ∈{ 1, …, K }


    Assume that {ξ̃^(t) = ξ̃_k^(t)} are independent sample points from the set {ξ_i, i = 1, …, I } taken with probabilities { p_i, i = 1, …, I }:
    
    ρ_t > 0,   ∑_t=0^∞ρ_t = ∞,   ∑_t=0^∞ρ_t^2 < ∞


    Let F(Y^*) denote the set of values of F on critical (stationary) points Y^* of problem (<ref>), where Y^* = { y = (y_1, …, y_K): ∂ F(y) ∈ N_Y (y_1) ×…× N_Y (y_K) } and N_Y (y_k) represents the normal cone to the set Y at point y_k. If F(Y^*) does not contain intervals and the sequence { y^(t)} is bounded, then { y^(t)} converges to a connected component of Y^*, and the sequence { F(y^(t)) } has a limit.




 §.§ Adaptive Stochastic Quantization
 
The minimization of the objective function (<ref>) is a non-smooth, non-convex, multiextremal, large-scale stochastic optimization problem. Although the parameter update recurrent sequence based on SGD (<ref>) can converge under conditions (<ref>), Qian et al. <cit.> demonstrated that the variance of gradient oscillations increases proportionally to the size of training samples:

    𝕍 (g_ℬ_k) ∝I^2/b𝕍 (g_k)

where 𝕍 represents the variance over a set, g_ℬ_k = 1/b∑_i=1^b g_i (ξ_ℬ_i) is the averaged gradient value over a subset ξ_ℬ_i⊂Ξ, and b = | ξ_ℬ_i |. These gradient oscillations reduce the algorithm's stability and slow down the convergence speed. While strategies such as manually tuned learning rate ρ > 0, annealing schedules <cit.>, or averaged gradient over a subset can improve convergence stability, the slow convergence speed in high-dimensional models <cit.> remains a significant drawback of the SGD algorithm.

Polyak <cit.> proposed the Momentum Gradient Descent (or the ”Heavy Ball Method”) as an alternative modification to the SGD by introducing an acceleration multiplier 0 < γ < 1 to the recurrent sequence (<ref>), using a physical analogy of the motion of a body under the force of friction:

    y_k+1 = y_k + γ (y_k - y_k-1) - ρ_k g_k^*


Nesterov <cit.> further improved the modified recurrent sequence (<ref>) by introducing an extrapolation step for parameter estimation (Nesterov Accelerated Gradient or NAG):

    ỹ_k = y_k - ρ_k g_k^*,    y_k+1 = ỹ_k + γ (ỹ_k - ỹ_k-1)


Although modifications (<ref>) and (<ref>) can improve convergence speed, they often encounter the vanishing gradient problem on sparse data <cit.>. The root cause is the fixed learning rate value, which performs equal updates for both significant and insignificant model parameters. Duchi et al. <cit.> address this issue by introducing an adaptive learning rate ρ̃_k = ρ_k / √(G_k + ε), where the hyperparameter value is normalized over the accumulated gradient value to increase the update for more significant parameters (AdaGrad):

    y_k+1 = y_k - ρ_k/√(G_k + ε) g_k^*

where G_k = G_k-1 + g_k^*^2 is a linear combination of accumulated gradients from previous iterations, and ε≪ 10^-8 is a denominator smoothing term. While approach (<ref>) solves the convergence issue on sparse data, it introduces the problem of uncontrollable vanishing of the learning rate with each iteration, i.e., lim_k →∞ | ρ̃_k | = 0. Tieleman et al. <cit.> proposed another approach (RMSProp) for accumulated gradient normalization using a moving average G_k = β G_k-1 + (1 - β) g_k^*^2, which substitutes the denominator G_k with a stochastic approximation of the expected value 𝔼 G_k to control learning rate vanishing with an averaging multiplier 0 < β < 1.

Kingma et al. <cit.> introduced a further modification to (<ref>) by adding adaptive estimation of the gradient value g_k^* (ADAM):

    m_k = β_1 m_k-1 + (1 - β_1) g_k 
    
            v_k = β_2 v_k-1 + (1 - β_2) g_k^2 
    
            y_k+1 = y_k - ρ_k/√(v_k + ε) m_k

where m_k is the adaptive first moment (expected value) estimation, v_k is the adaptive second moment (variance) estimation, and 0 < β_1 < 1, 0 < β_2 < 1 are averaging multipliers. It is important to note that the values m_i and v_i may be biased (i.e., the expected value of the parameter does not equal the value itself), which can cause unexpected behavior in the oscillation's variance. The authors proposed corrected estimations for (<ref>) as:

    m̅_k = m_k/1 - β_1,   v̅_k = v_k/1 - β_2


Norkin et al. <cit.> provide an overview of these adaptive parameter update strategies and present a detailed comparison of their convergence speed in various problem settings.

    

§ TRADITIONAL CLUSTERING MODEL


The optimal clustering problem seeks to find K cluster centers { y_1, ..., y_K } that minimize the sum of distances from each point ξ to the nearest center. Lloyd's K-Means algorithm <cit.> is a prominent method for solving quantization and clustering problems, with numerous extensions <cit.>. Bottou and Bengio <cit.> interpreted the K-Means algorithm as an analogue of Newton's method and proposed several stochastic gradient descent algorithms for optimal clustering. These stochastic K-Means algorithms use only a subset of possible ξ values or even a single element at each iteration.

     <cit.>. K-Means iterative algorithm starts with the set of current cluster centers { y_k^t,  k = 1, ..., K }, the feature set {ξ_i,  i = 1, ..., I } is subdivided into K non-intersecting groups { I_1^t, ..., I_K^t }: point ξ_i belongs to group I_s^t if
    
    ξ_i - y_s^t ^2 = min_1 ≤ k ≤ K{ξ_i - y_k^t ^2 },
    						    r>0.

    Denote N_k^t the number of points in group I_k^t, N_k^t = 0 if I_k^t = ∅, ∑_k=1^K N_k^t = I. Remark that I_k^t and N_k^t depend on y^t. K-Means iteratively evaluates next cluster centers y^t+1 with the estimation:
    
    y_k^t + 1 = 1/N_k^t∑_i ∈ I_k^tξ_i,  k = 1, ..., K;  t = 0, 1, ...

    We can represent these vectors as
    
    y_k^t + 1   =    y_k^t - 1/N_k^t∑_i ∈ I_k^t (y_k^t - ξ_i) = y_k^t - 1/N_k^t∑_i ∈ I_k^t y_k^t + 1/N_k^t∑_i ∈ I_k^tξ_i 
       =   1/N_k^t∑_i ∈ I_k^tξ_i,  k = 1, ..., K;  t = 0, 1, ...



In <cit.> the form of K-Means algorithm (<ref>) was connected to Newton's step for solving at each iteration the smooth quadratic problem

    min_y_1,…,y_K[F^t(y) = 1/2∑_k=1^K ∑_i ∈ I_k^t || ξ_i - y_k ||^2].

with block diagonal Hessian and with diagonal number 1/N_k^t in block k. Moreover, it is easy to see that (<ref>) is the exact analytical solution of the unconstrained quadratic minimization problem (<ref>) under fixed partition {I_1^t,…,I_K^t} of the index set {1,…,I}. In that paper it was also considered stochastic batch and online version of the stochastic gradient methods with learning rate 1/t + 1 for solving a sequence of problems (<ref>) but without rigorous convergence analysis.

The initial positions of cluster centers { y_k^0,  k = 1, ..., K } are set either at random among {ξ_i,  i = 1, ..., I } or according the algorithm K-Means++ <cit.>. With K-Means++ initialization strategy, the rate of convergence to local optimum is estimated to be 𝔼 [F] ≤ 8(ln k + 2 ) F^*, where F^* - an optimal solution <cit.>. Assume { y_1^0, ..., y_k^0 } (k<K) initial cluster centers have already been chosen. The next center y_k+1^0  (k+1<K) is sampled from the set {ξ_i,  i = 1, ..., I } with probabilities:

    q_j = min_1 ≤ s ≤ k || ξ_j - y_s^0 ||^2/∑_i=1^I min_1 ≤ s ≤ k || ξ_i - y_s^0 ||^2


The next positions of the cluster centers { y_1^t, ..., y_K^t } for t > 0 are calculated as in the original Lloyd algorithm <cit.>, the expectation-maximization (EM) approach to K-Means algorithm. Consider problem for r = 2: 

    F(y) = F(y_1, ..., y_K) = ∑_i=1^I p_i min_1 ≤ k ≤ K || y_k - ξ_i ||^2 = 𝔼_i ∼ pmin_1 ≤ k ≤ K || y_k - ξ_i ||^2


Given objective function (<ref>) is a generalized differentiable function <cit.>, and to find its optima we utilize a generalized gradient set is calculated by the chain rule:

    ∂ F(y) = ∑_i=1^I p_i ·conv_k ∈ K_i{ 2 (y_k - ξ_i) },  K_i = _1 ≤ k ≤ K || y_k - ξ_i ||


And its some generalized gradient is the compound vector g(y) = (g_k(y), k = 1, ..., K):

    g_k(y) = ∑_i=1^I 2 p_i (y_k_i - ξ_i),  k_i ∈ K_i


Sculley <cit.> addresses the limitation of Lloyd algorithm with update rule (<ref>), highlighting that objective function F(y) calculation is expensive for large datasets, due to O(K · I · d) time complexity for a given feature set {ξ_i }, where d - dimensionality of each sample ξ_i. The author proposed a solution by introducing a Mini-Batch K-Means modification to the Lloyd algorithm, where the set of points {ξ_i, i = 1, ..., I } is subdivided into non-intersecting subsets I_k, k = 1, ..., K,  such that i ∈ I_k if y_k - ξ_i  = min_k^'∈{ 1, ..., K } y_k^' - ξ_i.  Some I_k may occur empty, for example, if 
    max_1 ≤ i ≤ Imin_1 ≤ k ≤ Ky_k - ξ_i < min_1 ≤ i ≤ Iy_k - ξ_i.
  In this modifications the generalized gradients of function F(y) is g(y) = (g_1(y), ..., g_K(y)), where:

    g_k(y) = ∑_i ∈ I_k 2 p_i (y_k_i - ξ_i),     I_k ≠∅
    
            0,     I_k = ∅


The standard generalized gradient method for solving problem (<ref>) takes on the form (for p_i = 1/I, k = 1, ..., K, t = 0, 1, ...):

    y_k^t+1 = y_k^t - ρ_t g_k(y^t) = 
            y_k^t - ρ_t 2/I∑_i ∈ I_k (y_k^t - ξ_i),     I_k ≠∅
    
            0,     I_k = ∅


Recent studies <cit.> have examined stochastic K-Means algorithms as methods for solving corresponding non-convex, non-smooth stochastic optimization problems. However, the convergence properties of these algorithms lack rigorous validation. Let N_k^t is the number of elements in I_k ≠∅ at iteration t. If we chose ρ_t dependent on k, namely, ρ_tk = 0.5 I/N_k^t, then process (<ref>) becomes identical to K-Means one (<ref>). Here ρ_tk≤ 0.5 I and can be rather large that does not guarantee convergence of (<ref>). A more general choice can be ρ_tk = ρ_t I/N_k^t with ρ_t ≥ 0 satisfying conditions  
(<ref>) and thus lim_t →∞max_k ρ_tk = 0.



 §.§ Modifications of K-Means algorithm


Robust clustering model assumes solving problem (<ref>), (<ref>) with parameter r<2. Such choice of parameter r makes the quantization and clustering model more robust to outliers. However, then stochastic generalized gradients of the objective function should be calculated by formula (<ref>) and the stochastic clustering algorithm takes form (<ref>).

One also can consider a complement to the sequence (<ref>) by the Cesàro trajectory averaging:

    y̅_k^t+1 = (1 - σ_t+1) y̅_k^t + σ_t+1 y_k^t+1, σ_t+1 = ρ_t+1/∑_s=1^t+1ρ_s,  k = 1, ..., K


Conditions of convergence for this averaged sequence were studied in <cit.>, in particular, they admit learning rate ρ_t proportional to 1 / √(t+1). A similar approach for K-Means generated sequences (<ref>) aims to average sequence by the feature set size N_k:

    ỹ_k^t+1 = (1 - σ̃_t+1) ỹ_k^t + σ̃_t+1 y_k^t+1, σ̃_k, t+1 = 1/N_k^t+1 / ∑_s=1^t+11/N_k^s,  k = 1, ..., K


The standard K-Means algorithm requires finding the nearest cluster center _1 ≤ k ≤ K y_k - ξ_i to each point {ξ_i,  i = 1, ..., I }. This can be a time consuming operation in case of very large number I. Moreover, points ξ_i may be sampled sequentially from some continuous distribution and thus the sample {ξ_i,  i = 1, 2, ... } can be potentially arbitrary large. The Stochastic Quantization algorithm (<ref>) uses only one sampled point ξ̃^t at each iteration and thus only one problem _1 ≤ k ≤ K y_k - ξ̃^t is solved at iteration t. But one can use a batch of m such points {ξ̃^t_i_1^t, ..., ξ̃^t_i_m^t} instead of the whole sample {ξ_i }, m < I. 

    

§ DATA ENCODING WITH TRIPLET NETWORK


Our proposed semi-supervised classification approach addresses a critical challenge in clustering high-dimensional data by integrating Contrastive Learning techniques <cit.> with Stochastic Quantization. While the Stochastic Quantization algorithm (<ref>) effectively mitigates scalability issues for large datasets, it remains susceptible to the ”curse of dimensionality,” a phenomenon common to clustering methods that rely on distance minimization (e.g., K-Means and K-Means++). Kriegel et al. <cit.> elucidated this phenomenon, demonstrating that the concept of distance loses its discriminative power in high-dimensional spaces. Specifically, the distinction between the nearest and farthest points becomes increasingly negligible:

    lim_n →∞max d(ξ_i, y_k) - min d(ξ_i, y_k)/min d(ξ_i, y_k) = 0.


Our study focuses on high-dimensional data in the form of partially labeled handwritten digit images <cit.>. However, it is important to note that this approach is not limited to image data and can be applied to other high-dimensional data types, such as text documents <cit.>. While efficient dimensionality reduction algorithms like Principal Component Analysis (PCA) <cit.> exist, they are primarily applicable to mapping data between two continuous spaces. In contrast, our objective necessitates an algorithm that learns similarity features from discrete datasets and projects them onto a metric space where similar elements are grouped into clusters, a process known as similarity learning.

Recent research <cit.> has employed a Triplet Network architecture to learn features from high-dimensional discrete image data and encode them into low-dimensional representations in the latent space. The authors proposed a semi-supervised learning approach where the Triplet Network is trained on a labeled subset of data to encode them into latent representations in ℝ^n, and subsequently used to project the remaining unlabeled fraction onto the same latent space. This approach significantly reduces the time and labor required for data annotation without compromising accuracy.

The Triplet Network, introduced by <cit.>, is a modification of the Contrastive Learning framework <cit.>. Its core idea is to train the model using triplets of samples:

    
  * An anchor sample ξ_i: a randomly sampled element from the feature set Ξ
    
  * A positive sample ξ^+_i: an element with a label similar to the anchor ξ_i
    
  * A negative sample ξ^-_i: an element with a label different from the anchor ξ_i


Unlike traditional Contrastive Learning, which compares only positive ξ^+_i and negative ξ^-_i samples, the Triplet Network learns to minimize the distance between the anchor and positive samples while maximizing the distance between the anchor and negative samples. This is achieved using the triplet loss objective function (see Fig. <ref>):

    ℒ_triplet(θ) = max (0, d(f_θ(ξ_i), f_θ(ξ^+_i)) - d(f_θ(ξ_i), f_θ(ξ^-_i)) + α)

where f_θ: Ξ→ℝ^n is a parameterized abstract operator mapping discrete elements Ξ into latent representations (in our case, a Triplet Network with weights θ), d: [ℝ^n, ℝ^n] →ℝ is a distance metric between samples, and α is a margin hyperparameter enforcing a minimum separation between positive and negative pairs. Analogous to the Stochastic Quantization distance metric (<ref>), we employed the Euclidean norm l_2 for d(ξ_i, ξ_j) in (<ref>).


In our research, we utilized a Convolutional Network architecture as f_θ, as proposed by <cit.>. The detailed overview of the architecture, its training using the Backpropagation algorithm, and accuracy evaluation are beyond the scope of this paper; <cit.>, <cit.>, and <cit.> provide extensive coverage of these topics. However, notice that function (<ref>) is non-smooth and non-convex and the standard Backpropagation technique is not validated for such case. The extension of this technique to the non-smooth non-convex case was made in <cit.>.

Regarding triplet mining strategies for (ξ_i, ξ^+_i, ξ^-_i), it is crucial to select an approach that produces the most informative gradients for the objective function (<ref>). Paper <cit.> discusses various online triplet mining strategies, which select triplets within each batch of a training set on each iteration. We employed the semi-hard triplet mining strategy, which chooses an anchor-negative pair that is farther than the anchor-positive pair but within the margin α:

    ξ^-_i = _ξ: C(ξ) ≠ C(ξ_i) 
     d(f_θ(ξ_i), f_θ(ξ)) > d(f_θ(ξ_i), f_θ(ξ^+_i)) d(f_θ(ξ_i), f_θ(ξ))

where C(ξ) denotes the label of an element ξ.

By applying ideas from <cit.>, we can utilize the encoded latent representations of the Triplet Network to train a Stochastic Quantization (<ref>) algorithm. This novel approach enables us to solve supervised or semi-supervised learning problems of classification on high-dimensional data. The semi-supervised learning process using the combined algorithm, assuming we have a labeled subset ξ⊂Ξ and remaining unlabeled data ξ̅⊂Ξ, ξ̅∩ξ = ∅, can be summarized as follows:

    
  * Train a Triplet Network f_θ on labeled data ξ and produce encoded latent representations space Ξ̅⊂ℝ^n
    
  * Utilize the trained Triplet Network f_θ to project the remaining unlabeled data ξ̅ onto the same latent representation space Ξ̅
    
  * Employ both labeled and unlabeled latent representations to train a Stochastic Quantization (<ref>) algorithm

    

§ NUMERICAL EXPERIMENTS


To implement and train the Triplet Network, we utilized PyTorch 2.0 <cit.>, a framework designed for high-performance parallel computations on accelerated hardware. The Stochastic Quantization algorithm was implemented using the high-level API of Scikit-learn <cit.>, ensuring compatibility with other package components (e.g., model training and evaluation), while leveraging NumPy <cit.> for efficient tensor computations on CPU. All figures presented in this study were generated using Matplotlib <cit.>. The source code and experimental results are publicly available in our GitHub repository <cit.>.

For our experiments, we employed the original MNIST handwritten digit dataset <cit.> (see Fig. <ref>), comprising 60,000 grayscale images of handwritten digits with a resolution of 28x28 pixels, each associated with a class label from 0 to 9. Additionally, the dataset includes a corresponding set of 10,000 test images with their respective labels. It is noteworthy that we did not apply any data augmentation or preprocessing techniques to either the training or test datasets.


We approached the image classification task as a semi-supervised learning problem, training models on varying fractions of labeled training data (10%, 30%, 60%, 80%, and 100%). The training dataset was split using uniform sampling into labeled and unlabeled portions according to the specified percentages. For the Triplet Network, we employed a Convolutional Neural Network architecture consisting of two Convolutional Layers with 3x3 filters (feature map dimensions of 32 and 64, respectively, with stride=1 and padding=1), followed by 2x2 Max-Pooling layers, and two Dense Layers. ReLU (non-differentiable) activation functions were used throughout the network. Together with non-smooth Triplet loss function this makes the learning problem highly non-smooth and non-convex (see discussion of this issue in <cit.>). We trained separate Triplet Network models for each labeled data fraction with the following hyperparameters: 50 epochs, batch size of 1000 ×fraction, learning rate ρ = 10^-3, and l_2 regularization rate λ = 10^-5. For the triplet loss (<ref>) and triplet mining (<ref>), we set the margin hyperparameter α = 1.0. To facilitate meaningful feature capture while enabling visualization, we chose a latent space dimensionality of ℝ^3 (see Fig. <ref>).


The Triplet Network was used to project latent representations onto ℝ^3 space from both labeled and unlabeled training data. These representations were then used to train the Stochastic Quantization algorithm as an unsupervised learning model. For each set of latent representations corresponding to different labeled data fractions, we trained the Stochastic Quantization algorithm and its adaptive variants (<ref>)–(<ref>) from subsection <ref>. We employed the K-Means++ initialization strategy (<ref>) for all variants, with a rank hyperparameter r = 3. To ensure convergence, we used different learning rates for each variant: ρ = 0.001 for SGD, Momentum, and NAG; ρ = 0.9 for AdaGrad; and ρ = 0.01 for RMSProp and ADAM. With these hyperparameters, all Stochastic Quantization variants converged to the global optima (see Fig. <ref>), with most converging on the first iteration (see Fig. <ref>).



The accuracy of the trained classification models, combining Triplet Network and Stochastic Quantization, was evaluated using the F1-score metric <cit.> for weighted multi-label classification. Our experiments demonstrated that our approach achieved results comparable to state-of-the-art performance with Triplet Network and Siamese Network, as reported in <cit.>, even with limited labeled data:

    

§ CONCLUSIONS


In this paper, we introduced a novel approach to solving semi-supervised learning problems by combining Contrastive Learning with the Stochastic Quantization algorithm. Our robust solution addresses the challenge of scalability in large datasets for clustering problems, while also mitigating the ”curse of dimensionality” phenomenon through the integration of the Triplet Network. Although we introduced modifications to the Stochastic Quantization algorithm by incorporating an adaptive learning rate, there is potential for further enhancement by developing an alternative modification using the finite-difference algorithm proposed in <cit.>, which will be the focus of future research. The semi-supervised nature of the problems addressed in this paper arises from the necessity of using labeled data to train the Triplet Network. By employing other Deep Neural Network architectures to produce low-dimensional representations in the latent space, the semi-supervised approach of Stochastic Quantization with Encoding could be extended to unsupervised learning. Additionally, alternative objective functions for contrastive learning, such as N-pair loss <cit.> or center loss <cit.>, offer further avenues for exploration in future studies. The results obtained in this paper demonstrate that the proposed approach enables researchers to train classification models on high-dimensional, partially labeled data, significantly reducing the time and labor required for data annotation without substantially compromising accuracy.


    splncs04
    
