{'fraction_non_alphanumeric': 0.045454545454545456, 'fraction_numerical': 0.017686496222544652, 'mean_word_length': 4.440255851932498, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Two-stage Classification Method for High-dimensional Data and Point Clouds

High-dimensional data classification is a fundamental task in machine learning and imaging science. In this paper, we propose a two-stage multiphase semi-supervised classification method for classifying high-dimensional data and unstructured point clouds. To begin with, a fuzzy classification method such as the standard support vector machine is used to generate a warm initialization. We then apply a two-stage approach named SaT (smoothing and thresholding) to improve the classification. In the first stage, an unconstraint convex variational model is implemented to purify and smooth the initialization, followed by the second stage which is to project the smoothed partition obtained at stage one to a binary partition. These two stages can be repeated, with the latest result as a new initialization, to keep improving the classification quality. We show that the convex model of the smoothing stage has a unique solution and can be solved by a specifically designed primal-dual algorithm whose convergence is guaranteed. We test our method and compare it with the state-of-the-art methods on several benchmark data sets. The experimental results demonstrate clearly that our method is superior in both the classification accuracy and computation speed for high-dimensional data and point clouds.

1. Introduction. Data sets classification is a fundamental task in remote sensing, machine learning, computer vision, and imaging science [1,56,43,52,38]. The task, simply speaking, is to group the given data into different classes such that, on one hand, data points within the same class shares similar characteristics (e.g. distance, edges, intensities, colors, and textures); on the other hand, pairs of different classes are as dissimilar as possible with respect to certain features. In this paper, we focus on the task of multi-class semi-supervised classification. The total number of classes K of the given data sets is assumed to be known, and a few samples, namely the training points, in each class have been labeled. The goal is therefore to infer the labels of the remaining data points using the knowledge of the labeled ones.
For data classification, previous methods are generally based on graphical models, see e.g. [1,47,56], and references therein. In a weighted graph, the data points are vertices and the edge weights signify the affinity or similarity between pairs of data points, where the larger the edge weight is, the closer or more similar the two vertices are. The basic assumption for data classification is that vertices in the graph that are connected by edges with large weight should belong to the same class. Since a fully connected graph is dense and has the size as large as the square of the number of vertices, it is computationally expensive to work on it directly. In order to circumvent this, some cleverly designed approximations have been developed. For example, in [33,44], spectral approaches are proposed to efficiently calculate the eigendecomposition of a dense graph Laplacian. In [32,42], the nearest neighbor strategy was adopted to build up a sparse graph where most of its entries are zero, and therefore it is computationally efficient.
In the literature, various studies for semi-supervised classification have been performed by computing a local minimizer of some non-convex energy functional or minimizing a relevant convex relaxation. To name just a few, we have the diffuse interface approaches using phase field representation based on partial differential equation techniques [4,40], the MBO scheme established for solving the diffusion equation [45,44,33], and the variational methods based on graph cut [1,56]. In particular, in [1], the convex relaxation models and special constraints on class sizes were investigated. In [56], some novelty region-force terms were introduced in the variational models to enforce the affinity between vertices and training samples. To the best of our knowledge, all these proposed variational models have the so-called no vacuum and overlap constraint on the labeling functions, which gives rise to non-convex model with NP-hard issues. By allowing labeling functions to take values in the unit simplex, the original NP-hard combinatorial problem is rephrased into a continuous setting, see e.g. [1,56,43,52,57,36,10,18] for various continuous relaxation techniques (e.g. the ones based on solving the eigenvalue problem, convex approximation, or non-linear optimization) and references therein.
Image segmentation can also be viewed as a special case of the data classification problem [52,8], since the pixels in an image can be treated as individual points. Various studies and many algorithms have been considered for image segmentation. In particular, variational methods are among the most successful image segmentation techniques, see e.g. [46,30,13,16,31,2,9]. The Mumford-Shah model [46], one of the most important variational segmentation models, was proposed to find piecewise smooth representations of different segments. It is, however, difficult to solve since the model is non-convex and non-smooth. Then substantially rich follow-up works were conducted, and many of them considered compromise techniques such as: (i) simplifying the original complex model, e.g. finding piecewise constant solutions instead of piecewise smooth solutions [27,41,55]); (ii) performing convex approximations, e.g. using convex regularization terms like total variation [50,22]; or (iii) using the smoothing and thresholding (SAT) segmentation methodology [18,17,16,21,25]; for more details refer to e.g. [26,35,11,39,26,48,49,58,5] and references therein. Moreover, various applications were put forward for instance in optical flow [19], tomographic imaging [3], and medical imaging [59,14,15,12,51,20].
In this paper, we propose a multi-class semi-supervised data classification method based on the SaT segmentation methodology [18,17,16,21,25]. It has been shown to be very promising in terms of segmentation quality and computation speed for images corrupted by many different types of blurs and noises. Briefly speaking, the SaT methodology includes two main steps: the first step is to obtain a smooth approximation of the given image through minimizing some convex models; and the second step is to get the segmentation results by thresholding the smooth approximation, e.g. using thresholds determined by the K-means algorithm [37]. Since the models used are convex, the non-convex and NP-hard issues in many existing variational segmentation methods (e.g. the Mumford-Shah model and piecewise constant Mumford-Shah model mentioned above) were naturally avoided.
Our proposed data classification method mainly contains two stages with a warm initialization. The warm initialization is a fuzzy classification result which can be generated by any standard classification methods such as the support vector machine (SVM) [29]; or by labeling the given data randomly if no proper method is available for the given data (e.g. the data set is too large). Its accuracy is not critical since our proposed method will improve the accuracy significantly from this starting point.
With the warm initialization, the first stage of our method is to find a set of smooth labeling functions, where each gives the probability of every point being in a particular class. They are obtained by minimizing a properly-chosen convex objective functional. In detail, the convex objective functional contains K independent convex sub-minimization problems, where each corresponds to one labeling function, with no constraints between these K labeling functions. For each sub-minimization problem, the model is formed by three terms: (i) the data fidelity term restricting the distance between the smooth labeling function and the initialization; (ii) the graph Laplacian ( 2 -norm) term, and (iii) the total variation ( 1 -norm) built on the graph of the given data. The graph Laplacian and the total variation terms regularize the labeling functions to be smooth but at the same time close to a representation on the unit simplex.
After obtaining the set of labeling functions, the second stage of our method is just to project the fuzzy classification results obtained at stage one onto the unit simplex to obtain a binary classification result. This step can be done straightforwardly. To improve the classification accuracy, these two stages can be repeated iteratively, where at each iteration the result at the previous iteration is used as a new initialization.
The main advantage of our proposed method is twofold. First, it performs outstandingly in computation speed, since the proposed model is convex and the K sub-minimization problems are independent with each other (with no constraint on the K labeling functions). The parallelism strategy can be applied straightforwardly to improve computation performance further. On the contrary, the standard start-ofthe-art variational data classification methods e.g. [56,33,40] have the constraint on unit simplex in their minimization models, so that the non-convex or NP-hard issues can affect seriously the efficiency of these methods, even though some convex relaxations may be applied. Secondly, our method is generally superior in classification accuracy, due to its flexibility of merging the warm initialization and the two-stage iterations which are tractable and manage to improve the accuracy gradually. Note again that we are solving a convex model in the first stage of each iteration, which guarantees a unique global minimizer. In contrast, there is however no guarantee that the results obtained by the standard start-of-the-art variational data classification methods e.g. [56,33,40] are global minimizers. The effectiveness of iterations in our proposed method will be shown in the experiments. For most cases, the clustering accuracy would be increased by a significant margin compared to the first initialization and generally outperforms the state-of-the-art variational classification methods.
The paper is organized as follows. In Section 2, we give the basic notation used throughout the paper. In Section 3, we present our method for data sets classification. In Section 4, we present the algorithm for solving the proposed model and its convergence proof. In Section 5, we test our method on benchmark data sets and compare it with the start-of-the-art methods. Conclusions are drawn in Section 6.
2. Basic notation. Let G = (V, E, w) be a weighted undirected graph representing a given point cloud, where V is the vertex set (in which each vertex represents a point) containing N vertices, E is the edge set consisting of pairs of vertices, and w : E → R + is the weight function defined on the edges in E. The weights w(x, y) on the edges (x, y) ∈ E measure the similarity between the two vertices x and y; the larger the weight is, the more similar (e.g. closer in distance) the pair of the vertices is.
There are many different ways to define the weight function. Let d(·, ·) be a distance metric. Several particularly popular definitions of weight functions are as follows: (i) radial basis function w(x, y) := exp(−d(x, y) 2 /(2ξ)), ∀(x, y) ∈ E, (2.1) for a prefixed constant ξ > 0; (ii) Zelnic-Manor and Perona weight function where var(·) denotes the local variance; and (iii) the cosine similarity w(x, y) := x, y x, x y, y , ∀(x, y) ∈ E, where ·, · is the inner product. Let W = (w(x, y)) (x,y)∈E ∈ R N ×N , the so-called affinity matrix, which is usually assumed to be a symmetric matrix with non-negative entries. Let D = (h(x, y)) (x,y)∈E ∈ R N ×N be the diagonal matrix, where its diagonal entries are equal to the sum of the entries on the same row in W , i.e. |w(x, y)(u(x) − u(y))|, (2.6) and the 2 -norm (also known as Dirichlet energy) Note, however, that working with the fully connected graph E-like the setting in (2.5), (2.6) and (2.7)-can be highly computational demanding. In order to reduce the computational burden, one often only considers the set of edges with large weights. In this paper, the k-nearest-neighbor (k-NN) of a point x, N (x), is used to replace the whole edge set starting from the point x in E. Besides the computational saving, one additional benefit of using k-NN graph is its capability to capture local property of points lying close to a manifold. With the k-NN graph, then the definitions in (2.5), (2.6) and (2.7) become (2.8) and respectively, see e.g. [40,56] for more detail.
3.1. Preliminary. Given a point cloud V containing N points in R M . We aim to partition V into K classes V 1 , · · · , V K based on their similarities (the points in the same class possess high similarity), with a set of training points In other words, we aim to assign the points in V \ T certain labels between 1 to K using the training set T in which the labels of points are known, and the partition satisfies no vacuum and overlap constraint: In the rest of the paper, we denote the points in V needed to be labeled as S = V \ T , and call S the test set in V . The constraint (3.1) can be described by a binary matrix function U : Clearly, the above definition yields is also known as the indicator constraint on the unit simplex. Since the binary representation in (3.2) generally requires solving a non-convex model with NP-hard issue, a common strategy-the convex unit simplex-is considered as an alternative Note, importantly, that the convex constraint (3.3) can overcome the NP-hard issue and make some subproblems convex, but generally the whole model can still be nonconvex. Therefore, solving a model with constraint (3.1), (3.2), or (3.3) can be time consuming, see e.g. [40,56] for more detail. If a result satisfying (3.3) is not completely binary, a common way to obtain an approximate binary solution satisfying (3.2) is to select the binary function as the nearest vertex in the unit simplex by the magnitude of the components, i.e.
Here, e i is the K-length unit normal vector which is 1 at the i-th component and 0 for all other components.

Proposed method.
In this section, we present our novel two-stage method for data (e.g. point cloud) classification based on the SaT strategy which has been validated very effective in image segmentation. Our method can be summarized briefly as follows: first, a classification result is obtained as a warm initialization by using a classical, fast, but need not be very accurate classification method such as SVM [29]; then, a proposed two-stage iteration scheme is implemented until no change in the labels of the test points could be made between consecutive iterations. Specifically, at the first stage, we propose to minimize a novel convex model free of constraint (like those in (3.1), (3.2) and (3.3)), to obtain a fuzzy partition, say U , while keeping the training labels unchanged; at the second stage, a binary result is obtained by just applying the binary rule in (3.4) directly on the fuzzy partition obtained at the first stage. This binary result could be the final classification result for the original classification problem, or, if necessary, be set as a new initialization to search a better one in the same manner. In the following, we give the details of each step.
Initialization. Given a point cloud V containing N points in R M and training set T containing N T points with correct labels, we use SVM, which is a standard and fast clustering method as an example, to obtain the first clustering. Let the partition matrix beÛ = (û 1 , · · · ,û K ) ∈ R N ×K , whereû j = (û j (x)) x∈V ∈ R N for j = 1, . . . , K. One could acquire an initialization by any other methods which have better performance than SVM. If no proper method is available (e.g. the data set is too large), then an initialization generated by setting labels to the test points randomly can be used as an alternative.
Stage one. Now we put forward our convex model to find a fuzzy partition U with initializationÛ = (û 1 , · · · ,û K ). It is where the first term is the data fidelity term constraining the fuzzy partition not far away from the initialization; the second term is related to ∇u 2 2 with graph Laplacian L; the last term is the total variation constructed on the graph; and α, β > 0 are regularization parameters. Specifically, the second term in (3.5) is used to impose smooth features on the labels of the points, and the last term is used to force the points with similar information to group together.
We emphasize that we already have the labels on the points in the training set T , withŪ = (ū 1 , · · · ,ū K ) ∈ R N T ×K being the partition matrix on T , whereū j = (ū j (x)) x∈T ∈ R N T for j = 1, . . . , K. Therefore, we only assign labels to points in the test set S, i.e. we haveû Letû Sj represent the part ofû j defined on the test set S, then we havê We use analogous notations for the partition matrix U = (u 1 , · · · , u K ), with In Section 4, (3.7) and (3.8) are going to be used to derive an efficient algorithm to solve (3.5). The following Theorem 3.1 proves that our proposed model (3.5) has a unique solution.
Proof. According to [7,Chapter 9], a strongly convex function has a unique minimum. The conclusion follows easily from the strong convexity of model (3.5).
Many algorithms can be used to solve model (3.5) efficiently due to the convexity of the model without constraint. For example, the split-Bregman algorithm [34], which is specifically devised for 1 regularized problems; the primal-dual algorithm [23], which is designed to solve general saddle point problems; and the powerful alternative, ADMM algorithm [6]. In particular, model (3.5) actually contains K independent sub-minimization problems, where each corresponds to a labeling function u j , and therefore the parallelism strategy is ideal to apply. This is an important advantage of our method for large data sets. The algorithm aspects to solve our proposed convex model (3.5) are detailed in Section 4.
Stage two. This stage is to project the fuzzy partition result U obtained at stage one to a binary partition. Here, formula (3.4) is applied to the fuzzy partition U to generate a binary partition, which naturally satisfies no vacuum and overlap constraint (3.1). We remark that compared to the computation time at stage one, the time at stage two is negligible.
Normally, the classification work is complete after we obtain a binary partition matrix at stage two. However, since the way of obtaining an initialization in our scheme is open, and the quality of the initialization could be poor, we suggest going back to stage one with the latest obtained partition as a new initialization and repeat the above two stages until no more change in the partition matrix is observed. More precisely, we set U asÛ and repeat Stages 1 and 2 again to obtain a new U . Then the final classification result is the converged stationary partition matrix, say U * . Moreover, to accelerate the convergence speed, we update β in (3.5) by a factor of 2 if we are to repeat the stages. This will obviously enforce the closeness between two consecutive clustering results during iterations, which will ensure the algorithm converges fast. We stop the algorithm when no changes are observed in the clustering result compared to the previous one. We remark that a few iterations (≈ 10) are generally enough in practice, see the experimental results in Section 5 for more detail.
Note, importantly, that our classification method here is totally different from other variational methods like [40,56] which need to minimize variational models with constraint like (3.1), (3.2), (3.3), or other kinds of constraints (e.g. minimum and maximum number of points imposed on individual classes V i ). Even though our proposed model (3.5) has no constraint, the final classification result of our method naturally satisfies no vacuum and overlap constraint (3.1). Therefore, our method is much easier to solve for each iteration. Our proposed method, namely SaT (smoothing and thresholding) method for high-dimensional data classification, is summarized in Algorithm 1. 4. Algorithm aspects. In this section, we present an algorithm to solve the proposed convex model (3.5) based on the primal-dual algorithm [23] which is briefly recalled below.

Algorithm 1 SaT method for high-dimensional data classification
4.1. Primal-dual algorithm. Let X i be a finite dimensional vector space equipped with a proper inner product ·, · Xi and norm · Xi , i = 1, 2. Let map K : X 1 → X 2 be a bounded linear operator. The primal-dual algorithm is, generally speaking, to solve the following saddle-point problem where G : X 1 → [0, +∞], F : X 2 → [0, +∞] are proper, convex, lower-semicontinuous functions, and F * represents the convex conjugate of F. Given proper initializations, the primal-dual algorithm to solve (4.1) can be summarized in the following iterative way of updating the primal and dual variables 2) where θ ∈ [0, 1], τ, σ > 0 are algorithm parameters.

4.2.
Algorithm to solve our proposed model. We first define some useful notations which will be used to present our algorithm.

Preliminary.
For ease of explanation, in the following, when we say (i, j) ∈ E, the i and j represent the i-th and j-th vertices in E, respectively. Let (4.5) The graph laplacian L = D − W ∈ R N ×N can be decomposed as where is a matrix with only four nonzero entries which locate at positions (i, i), (i, j), (j, i) and (j, j).
Then the decomposition L in (4.6) can be rewritten as Note that, the terms (i,j)∈E a L ij and (i,j)∈E b L ij only have nonzero entries which are associated to the test set S and the training set T , respectively. Let According to (2.8), the gradient operator ∇ can be regarded as a linear transformation between R N and R N ×(k−1) (where k = |N (x)|). For a vector u j = (u Sj ,ū j ) defined in (3.8), let is an operator corresponding to the test set S, and H j is the gradient matrix corresponding to the training set T which is fixed sinceū j is fixed. Then, we have (4.16) Note, obviously, that solving the above model (4.16) is equivalent to solving K subminimization problems corresponding to each u Sj , j = 1, . . . , K, which means that our proposed model inherently benefits from the parallelism computation.
Using the definition of the 1 -norm given in (2.9), the conjugate of F j , F * j , can then be calculated as where P = {p ∈ R N ×(k−1) : p ∞ ≤ 1}, and χ P (p) is the characteristic function of set P with value 0 if p ∈ P , otherwise +∞. Using the primal-dual formulation in (4.1) with the definitions of G j and F * j respectively given in (4.17) and (4.19), then the minimization problem (4.16) corresponding to each u Sj can be reformulated as To apply the primal-dual method, it remains to compute (I+σ∂F * j ) −1 and (I + τ ∂G j ) −1 . Firstly, for ∀x ∈ R N ×(k−1) , we have where the operator ι P (·) is the pointwise projection operator onto the set P , i.e., ∀p ∈ R, Secondly, for ∀x ∈ R N −N T , we have Using the definition of G j (u Sj ) given in (4.17), problem (4.23) becomes solving the following linear system Since (αL + βI + 1 τ I) is positive definite, the above linear system can be solved efficiently by e.g. conjugate gradient method [24].
Finally, by exploiting the strong convexity of G j , ∀1 ≤ j ≤ K, which is shown in the next lemma, [23] suggests that we could adaptively modify σ, τ to accelerate the convergence or the primal-dual method.
Lemma 4.1. The functions G j , ∀1 ≤ j ≤ K are strongly convex with parameter β.
Proof. For simplicity, we omit the subscript j and S j in the following proof. First, by (4.12), L S is semi-positive definite. Therefore, ( α 2 u L S u + αu L 3ū ) is convex. Now the strong convexity of G follows from the fact that the remaining term in (4.17), which is β 2 u −û 2 2 , is strongly convex with parameter β. The algorithm solving our proposed classification model (4.16) (i.e. model (3.5)) is summarized in Algorithm 2, and its convergence proof is given in Theorem 4.2 below. For each sub-minimization problem, the relative error between two consecutive iterations and/or a given maximum iteration number can be used as stopping criteria to terminate the algorithm. Finally, we emphasize again that our method is quite suitable for parallelism since the K sub-minimization problems are independent with each other and therefore can be computed in parallel.

Now, we have
Therefore, we conclude that, the algorithm converges as long as we choose τ (0) , σ (0) > 0, such that τ (0) σ (0) < 1 N 2 (k−1) . 5. Numerical results. In this section, we evaluate the performance of our proposed method on four benchmark data sets-including Three Moon, COIL, Opt-Digits and MINST-for semi-supervised learning. Three Moon is a synthetic data set which has been used frequently e.g. [33,40,56]. The COIL, Opt-Digits, and MNIST data set can be found from the supplementary material of [28], the UCI machine learning repository 1 , and the MNIST Database of Handwritten Digits 2 , respectively. The basic properties of these test data sets are shown in Table 5.1. To implement our method, k-NN graphs are constructed for the test data sets, using the randomized kd-tree [53] to find the nearest neighbors with Euclidean distance as the metric. The radial basis function (2.1) is used to compute the weight matrix W , except for the MNIST data set where the Zelnic-Manor and Perona weight function (2.2) is used with eight closed neighbors. The training samples T -samples with labels known-are selected randomly from each test data set. The classification accuracy is defined as the percentage of correctly labeled data points.
For our proposed method, the regularization parameter β is fixed to 10 −4 for MNIST, 10 −5 for COIL, and 10 −2 for Three Moon, Opt-Digits. In practice, one could choose the value of β based on the accuracy of initialization. The better the initialized accuracy, the larger β we could choose. The regularization parameter α is set to 1 for Three Moon, Opt-Digits, 0.4 for MNIST, and 10 −2 for COIL. The accuracy of the proposed method can be improved further after fine-tuning the values of α and β for individual test data sets. All the codes were run on a MacBook with 2.8 GHz processor and 16 GB RAM, and Matlab 2017a.

Methods comparison.
As mentioned in previous sections, we use SVM method [29] to generate initializations for our proposed method. If it is not proper for a data set (e.g. very slow due to the large size of the data set), we could just use an initialization generated by assigning clustering labels randomly.
The SVM is a technique aiming to find the best hyperplane that separates data points of one class from the others. In practice, data may not be separable by a hyperplane. In that case, soft margin is used so that the hyperplane would separate many data points if not all. It is also common to kernelize data points, and then find separating hyperplane in the transformed space. The SVM method used in our experiments is trained with linear kernel.
We compare our proposed method with the state-of-the-art methods proposed recently, e.g. CVM [1], GL [33], MBO [33], TVRF [56], LapRF [56], LapRLS [54], MP [54], and SQ-Loss-I [28]. The code TVRF was provided by the authors and the parameters used in it were chosen by trial and error to give the best results. The classification accuracies of methods GL, MBO, LapRF, LapRLS, MP and SQ-Loss-I were taken from [1,56], in which methods CVM and TVRF were shown to be superior in most of the cases.

5.2.
Three moon data. The synthetic Three Moon data used here is constructed by following the way performed in [1,56] exactly. We briefly repeat the procedure as follows. First, generate three half circles in R 2 -two half top unit circles and one half bottom circle with radius 1.5 which are centered at (0, 0), (3, 0) and (1.5, 0.4), respectively. Then 500 points are uniformly sampled from each half circle and embedded into R 100 by appending zeros to the remaining dimensions. Finally, an i.i.d Gaussian noise with standard deviation 0.14 is added to each dimension of the data. An illustration of the first two dimensions of the Three Moon data is shown in Fig. 5.1 (a) where different colors are applied on each half circle. This is a three-class classification problem with the goal of classifying each half circle using a small number of supervised points from each class. This classification problem is challenging due to the noise and the high dimensionality of all the points with high similarity in R 98 . A k-NN graph with k = 10 is built for this data set, parameter σ = 3 is used in the Gaussian weight function, and the distance metric chosen is Euclidean metric for R 100 . We first test the methods using uniformly distributed supervised points, where a total number of 75 points are sampled uniformly from this data set as training points.
The accuracies of method TVRF and ours are obtained by running the methods ten times with randomly selected labeled samples, and taking the average of the accuracies. The accuracies of method CVM are copied from the original paper [1]. The accuracy comparison is reported in Table 5.2, which shows that our proposed method gives the highest accuracy; also, see Fig. 5.1 for visual validation of the results between methods of TVRF and ours. The average number of iterations taken for our proposed method is 3.8. Fig. 5.4 (a) gives the convergence history and partition  accuracy of our proposed method corresponding to iteration steps, which clearly shows the accuracy increment during iterations (note that the accuracy at iteration 0 is the result of the initialization which is obtained by SVM method). Table 5.7 reports the comparison in terms of computation time, which indicates the superior performance of our proposed method in computation speed. In the following, as a showcase, we test the methods using non-uniformly distributed supervised points, which is used to investigate the robustness of these methods on training points. In this case for the 75 training points, as an example, we respectively pick 5 points from the left and the bottom half circles, and pick the rest 65 points from the right half circle. This sampling is illustrated in Fig 5.2. The accuracies of TVRF and our method are shown in Table 5.3, from which we see clearly that the proposed method gives much higher accuracy. In particular, compared to the results in Table 5.2 using training points chosen uniformly, while the accuracy of TVRF method decreases by 0.8 percent, we observe only a very small decrease in our proposed method. This shows the robustness of our method with respect to the way that training points are selected. Note that in the case of training points chosen non-uniformly, the initialization obtained by SVM is poor, because of which more iterations are needed to converge for our method-average 12.0 iterations in 10 trials versus 3.3 iterations needed for the case of training points selected uniformly.

COIL data.
The benchmark COIL data comes from the Columbia object image library. It contains a set of color images of 100 different objects. These images, with size 128 × 128 each, are taken from different angles in steps of 5 degrees, i.e., 72 (= 360/5) images for each object. In the following, without loss of generality, we also call an image a point for ease of reference. The test data set here is constructed the same way as depicted in e.g. [1,56] and is briefly described as follows. First, the red channel of each image is down-sampled to 16 × 16 pixels by averaging over blocks of 8 × 8 pixels. Then, 24 out of the 100 objects are randomly selected, which amounts to 1728 (= 24 × 360/5) images. After that, these 24 objects are partitioned into six classes with four objects-288 images (= 4 × 72)-in each class. Finally, after discarding 38 images randomly from each class, a data set of 1500 images where 250 images in each of the six classes are constructed. To construct a graph, each image, which is a vector with length 241, is treated as a node on the graph, For accuracy test, a k-NN graph with k = 4 is built for this data set, parameter σ = 250 is used in the Gaussian weight function, and the distance metric chosen is Euclidean metric for R 241 . The training points, amount to 10% of the points, are selected randomly from the data set. Again, we run the test methods 10 times and compare the average accuracy. The resulting accuracies are listed in Table 5.4, showing that our method outperforms other methods. Moreover, the average number of iterations of our method is 12.2. Fig. 5.4 (b) gives the convergence history of our proposed method in partition accuracy corresponding to iterations, which again shows an increasing trend in accuracy. 5.4. MNIST data. The MNIST data set consists of 70,000 images of handwritten digits 0-9, where each image has a size of 28 × 28. Fig. 5.3 shows some images of the ten digits from the data set. Each image is a node on a constructed graph. The objective is to classify the data set into 10 disjoint classes corresponding to different digits. For accuracy test, a k-NN graph with k = 8 is built for this data set, and Zelnik-Manor and Perona weight function (2.2) is used to compute the weight matrix. The training 2500 (3.57%) points (images) are selected randomly from the total 70,000 points.
The experimental results of the test methods are obtained by running them 10 times with randomly selected training set with a fixed number of points 2500, and the average accuracy is computed for comparison. The accuracies of the test results are shown in Table 5.5, which indicates that our method is comparable to or better than the state-of-the-art methods compared here. Table 5.7 shows the computation time comparison, from which we again see that our method is very competitive in computation speed. The convergence history of our proposed method in partition accuracy corresponding to iterations is given in Fig. 5.4 (c), which also demonstrates a clear increasing trend in accuracy.

5.5.
Opt-Digits data. The Opt-Digits data set is constructed as follows. It contains 5620 bitmaps of handwritten digits (i.e. 0-9). Each bitmap has the size of 32 × 32 and is divided into non-overlapping blocks of 4 × 4, and then the number of "on" pixels are counted in each block. Therefore, each bitmap corresponds to a matrix of 8 × 8 where each element is an integer in [0, 16]. The classification problem is to partition the data set into 10 classes.
For accuracy test, a k-NN graph with k = 8 is built for this data set, parameter σ = 30 is used in the Gaussian weight function, and the distance metric chosen is Euclidean metric for R 64 . For the experiments on this data set, we generate three training sets respectively with the number of points 50, 100 and 150, which are all selected randomly. All the methods are run 10 times for each training set and the average accuracy is used for comparison. The quantitative results in accuracy are listed in Table 5.6, from which we see that our proposed method is consistently better than the state-of-the-art methods compared for all the cases. We also observe the improvement of the accuracy of these methods w.r.t. the increasing number of points in the training set. Finally, we show the convergence history of our proposed method in partition accuracy corresponding to iterations using the training set with 150 points in Fig. 5.4 (d), which again clearly shows an increasing trend in accuracy.   our proposed method outperforms the state-of-the-art methods for high-dimensional data and point clouds classification. Compared to the start-of-the-art variational classification models proposed e.g. in [1,56], in addition to the data fidelity term and 1 term (e.g. TV), our proposed model in (3.5) contains an additional 2 term on the labeling functions which is used to smooth the classification results so as to reduce the non-smooth artifact (the socalled staircase artifact in images) introduced by the 1 term. This is one reason that our method can generally achieve better results. Moreover, the warm initialization used in our method can also play a role to improve the classification quality. Apart from generating the initialization manually, any classification methods can practically be used to generate the initialization. Starting from the initialization, our proposed method can then be applied to achieve a better classification result by improving the accuracy iteratively. Theoretically speaking, the poorer the quality of the initialization, the more iterations are needed for our method. Nevertheless, we found that even for poor initializations (e.g. the ones generated randomly), 20 iterations are already enough to achieve competitive results. Generally, no more than 15 iterations are needed when using an initialization computed by standard classification methods (e.g. SVM).
Another distinction of our proposed model compared to the variational classification models in e.g. [1,56] is that there are no constraints on these labeling functions in our objective functional. In other words, in each iteration, we just need to find the minimizer of the objective functional corresponding to each labeling function, but these minimizers do not need to satisfy the constraint that their summation equal to 1. Therefore, the computation speed for every single iteration is improved in our method compared to other methods which have constraints. We emphasize again that, since minimizing each sub-problem with respect to each labeling function is irrelevant to minimizing the sub-problems with respect to other labeling functions, parallelism techniques can be used straightforwardly to further improve the computation performance of our algorithm; theoretically, we just require 1/K of the computation time needed for the non-parallelism scheme. This will be extremely important for large data sets. From Table 5.7, we see that, for all the computation time of our method, when considering the effect of parallel computing, our method should be able to outperform the state-of-the-art methods by a large margin.
6. Conclusions. In this paper, a two-stage multiphase semi-supervised method is proposed for classifying high-dimensional data or unstructured point clouds. The method is based on the SaT strategy which has been shown very effective for segmentation problems such as gray or color images corrupted by different degradations. Starting with a proper initialization which can be obtained by using any standard classification algorithm (e.g. SVM) or constructed by users, the first stage of our method is to solve a convex variational model without constraint. Most importantly, our proposed model is a lot easier to solve than the state-of-the-art variational models proposed recently (e.g. [1,56]) for point clouds classification problem since they all need no vacuum and overlap constraint (3.1) on the labeling functions in the unit simplex which could make their models to be non-convex. The second stage of our method is to find a binary partition via thresholding the smoothed result obtained from stage one. We proved that our proposed model has a unique solution and the derived primal-dual algorithm converges.
We tested our proposed method on four benchmark data sets and compared with the state-of-the-art methods. We also investigated the influence of the training sets selected uniformly and non-uniformly. For our method, different ways of generating initializations were implemented and validated. On the whole, the experimental results demonstrated that our method is superior in terms of classification accuracy and when parallel computing is considered, computation speed too. Therefore our method is an efficient and effective classification method for data sets like high-dimensional data or unstructured point clouds.