1 Introduction
In machine learning it is common to interpret each data point as a vector in Euclidean space [3]
. Such a discretisation is chosen because it allows for easy closed form solutions and fast computation, even with large datasets. However these methods ignore the fact that the data may not naturally fit into this assumption. In fact much of the data collected for practical machine learning are actually functions i.e. curves. For example financial data such as stock or commodity prices are functions of monetary value over time. Functional data have become increasingly important in many scientific and engineering research areas such as ECG (electrocardiogram) or EEG (Electroencephalography) in healthcare, biology data analysis, weather or climate data and motion trajectories from computer vision.
Analyzing functional data has been an emerging topic in statistical research [7, 13, 20, 21] and has attracted great attention from machine learning community in recent years [2, 15]
. One of important challenges in analyzing functional data for machine learning is to efficiently cluster and to learn better representations for functional data. Theoretically the underlying process for functional data is of infinite dimension, thus it is difficult to work with them with only finite samples available. A desired model for functional data is expected to properly and parsimoniously characterize the nature and variability hidden in the data. The classic functional principal component analysis (fPCA)
[17] is one of such examples to discover dominant modes of variation in the data. However fPCA may fail to capture patterns if the functional data are not well aligned in its domain. For time series, a special type of functional data, dynamic time warping (DTW) has long been proposed to compare time series based on shape and distortions (e.g., shifting and stretching) along the temporal axis [16, 24].Another important type of functional data is shape [23, 20]. Shape is an important characterizing feature for objects and in computer vision shape has been widely used for the purpose of object detection, tracking, classification, and recognition. In fact, a natural and popular representation for shape analysis is to parametrize boundaries of planar objects as 2D curves. In object recognition, images of the same object should be similar regardless of resolution, lighting, or orientation. Hence an efficient shape representation or shape analysis scheme must be invariant to scale, translation and rotation. A very useful shape representation is the squareroot velocity function (SRVF) representation [9, 20]. In general, the resulting SRVF of a continuous shape is square integrable, a welldefined Hilbert space where appropriate measurement can be applied, refer to Section 2 for more details. By acknowledging the true nature of the data we can develop more powerful methods that exploit features that would otherwise be ignored or lead to erroneous results with simple linear models.
Our intention in this study is to consider functional data clustering by accounting for the possible invariance in scaling/stretching, translation and rotation of functional data to help maintain shape characteristics. The focus of this paper is upon functional data where data sets consist of continuous real curves including shapes in Euclidean spaces. More specifically we propose a method of subspace analysis for functional data based on the idea developed in recent subspace clustering. The idea is to apply a feature mapping such as the aforementioned SRVF to the curves so that they are transformed onto the curve manifold, where the subspace analysis can be conducted based on the geometry on the manifold. In particular, we adapt the well known lowrank representation (LRR) framework [12] to deal with data that lie on the manifold of open curves by implementing the classical LRR in tangent spaces of the manifold [8, 25, 29].
LRR on Euclidean spaces [12] is closely related to several stateoftheart subspace analysis approaches such as Sparse Subspace Clustering (SSC) [6], Robust PCA (RPCA) [5] and lowrank Matrix Completion (MC) [28] methods. This mixture of subspaces model has naturally led to the development of subspace segmentation methods. Such methods aim to segment the data into clusters with each cluster corresponding to a unique subspace. More formally, given a data matrix of observed columnwise data samples , the objective of subspace clustering is to assign each data sample to its underlying subspace. The basic assumption is that the data within is drawn from a union of subspaces of dimensions .
The core of both SSC and LRR is to learn an affinity matrix for the given dataset and the learned affinity matrix will be pipelined to a spectral clustering method like nCUT
[19] to obtain the final subspace labels. To learn the affinity matrix, SSC relies on the self expressive property [6], which is thateach data point in a union of subspaces can be efficiently reconstructed by a linear combination of other points in the data.
In other words, each point can be written as a linear combination of the other points i.e. , where is a matrix of coefficients. Most methods however assume the data generation model , where is the observed data and is noise. Since it is difficult to separate the noise from the data the solution is to relax the selfexpressive model to , where is a fitting error and is different from .
Similarly LRR [12] exploits the self expressive property but attempts to learn the global subspace structure by computing the lowestrank representation of the set of data points. In other words, data points belonging to the same subspace should have similar coefficient patterns. In the presence of noise LRR attempts to minimise the following objective
(1) 
However rank minimisation is an intractable problem. Therefore LRR actually uses the nuclear norm
(sum of the matrix’s singular values) as the closest convex relation
(2) 
where is a placeholder for the norm most appropriate to the expected noise type. For example in the case of Gaussian noise the best choice is the norm i.e. and for sparse noise the norm should be used.
Both SSC and LRR rely on the linear self expressive property. This property is no longer available in the nonlinear manifold, e.g. the manifold of open curves as mentioned previously. To generalize LRR or SSC for data in the manifold space, we explicitly explore the underlying nonlinear data structure and utilize the techniques of exponential and logarithm mappings to bring data to a local linear space.
The rest of the paper is organized as follows. In Section 2, we review the preliminaries about the manifold of open curves and introduce the curve LowRank Representation (cLRR) model. Section 3 is dedicated to explaining an efficient algorithm for solving the optimization proposed in cLRR based on the linearized alternative direction method with adaptive penalty (LADMAP) and the algorithm convergence and complexity are also analyzed. In Section 4, the proposed model is assessed on both synthetic and real world databases against several stateoftheart methods. Finally, conclusions are discussed in Section 5.
2 LRR over the Curve Manifold
As previously discussed LRR is limited to a linear model and its current version can only be applied to vector data from a Euclidean space. Matrix in (1) or (2) encodes the affinity/similarity between data points. However this assumption is often unnatural and quite limiting. Much of the data encountered in real world is functional. In other words it exhibits a curve like structure over a domain. Euclidean linear models are unable to capture the nonlinear invariance embedded in each data point. For example in thermal infrared data of geological substances a curve may contain a key identifying feature such as a dip near a particular frequency. This dip may shift or vary position over time even for the same substance due to impurities. Under a linear vector model this variation may cause the vector to drastically move in the ambient Euclidean space and cause poor results. Or in other cases the feature may be elongated, shrunk or be subject to some nonuniformly warping or scaling. In all these cases the linear model will fail to accurately represent the nonlinear affinity in the data.
Exploring these unique nonlinear invariance in functional data is the focus of this paper. We now discuss how to adapt LRR (similar approach appliable to SSC) such that it easily accepts curve data and nonlinear relationships within clusters can be easily discovered.
2.1 The Curve Manifold
Given a smooth parameterized dimension curve , we represent it using he squareroot velocity function (SRVF) representation [9, 20], which is given by
The SRVF mapping transforms the original curve into a gradient based representation, which facilitates the comparing of the shape information.
In this paper, we focus on the set of open curves, e.g. the curves do not form a loop (). For handling general curves, we refer readers to [20]. The SRVF facilitates a measure and geometry bearing invariance to scaling, shifting and reparameterization in the curves domain. For example, all the translated curves from a curve will have the same SRVF. Robinson [18] proved that if the curve is absolutely continuous, then its SRVF is squareintegrable, i.e., is in a functional Hilbert space . Conversely for each , there exists a curve whose SRVF corresponds to . Thus the set is a welldefined representation space of all the curves. The most important advantage offered by the SRVF framework is that the natural and widely used measure on is invariant to the reparameterization. That is, for any two SRVFs and and a randomly chosen reparametrization function (nondecreasing) , we have
This property has been exploited in [2] for functional data clustering under the subspace clustering framework. Different from the work proposed in [2], we will adopt the newly developed LRR on manifolds framework to the model of curves LRR, see [8, 25, 29]. To see this, we introduce some more notation. Let be the set of all diffeomorphisms from to . This set collects all the reparametrization mappings. is a Lie group with the composition as the group operation and the identity mapping as the identity element. Then all the orbits together define the quotient manifold .
Without loss of generality, all curves are normalized to have unit length, i.e., . The SRVFs associated with these curves are elements of a unit hypersphere in the Hilbert space , i.e., . Therefore, under the curve normalization assumption, instead of , we consider the following unit hypersphere manifold
The manifold has some nice properties, see [1]. For any two points and in , a geodesic connecting them is given by ,
(3) 
where is the length of the geodesic. If we take derivative of w.r.t to , the tangent vector at is
(4) 
The above formula is regarded as the L͡ogarithm mapping on the manifold .
As we are concerned with the shape invariance, i.e., we need to additionally remove the shapepreserving transformations: rotation and curve reparametrization. The manifold concerning us is the quotient space of the manifold , where is the rotation group. Each element is an equivalent class defined by
Given any two points and in , a tangent representative [1] in the tangent space can be calculated in the following way, as suggested in [31, 22] based on (4),
(5) 
where is the representative of given by the welldefined algorithm in [31, 22] and . In fact, is the lifting representation of abstract tangent vector on at .
2.2 The Proposed Curve LRR
Given a set of unitlength curves , denote their SRVFs by such that and is a representative of the equivalent class . We cannot apply the standard LRR model (2) directly on the quotient manifold . This is because (2) indeed relies on the following individual linear combination
(6) 
which is invalid for ’s on . Note that can be explained as the affinity or similarity between data points and .
On any manifold, the tangent space at a given point is linearly local approximation to the manifold around the point and the linear combination is valid in the tangent space. This prompts us to replace the affinity relation in (6) by the following
(7) 
with the constraint to maintain consistency at different locations. The meaning of in (7) is the similarity between curves and via the “affinity” between tangent vectors and at the first order approximation accuracy. Each can be calculated by (5) and it is obvious that for any .
With all the ingredients at hand, we are fully equipped to propose the curve LRR (cLRR) model as follows
(8) 
where is the metric defined on the manifold, which is defined by the classic Hilbert metric on the tangent space.
Denote the th row of matrix and define
(9) 
Then with some algebraic manipulation we can rewrite the model (8) into the following simplified form,
(10) 
where .
Effectively this objective allows for similarity between curves to be measured in their tangent spaces. Our highly accurate segmentation results in Section 4 have demonstrated that this is an effective way to learn nonlinear similarity.
3 Optimisation
3.1 Algorithm
To solve the cLRR objective we use the Linearized Alternative Direction Method with Adaptive Penalty (LADMAP) [10, 11] . First take the Augmented Lagrangian of the objective (10)
(11)  
where is the Lagrangian multiplier (vector) corresponding to the equality constraint , is the matrix Frobebiusnorm, and we will update as well in the iterative algorithm to be introduced.
Denote by the function defined by (11) except for the first term . To solve (11), we adopt a linearization of at the current location in the iteration process, that is, we approximate by the following linearization with a proximal term
where is an approximate constant with a suggested value given by , and is a gradient matrix of at . Denote by
the 3order tensor whose
th front slice is given by . Let us define the matrix whose row is given by , then it is easy to show(12) 
Then (11) can be approximated by linearization and will be updated by the following
(13)  
Problem (13) admits a closed form solution by using SVD thresholding operator [4], given by
(14) 
where is the SVD of and is the Singular Value Thresholding (SVT) [4, 14] operator defined by
(15) 
The updating rule for
(16) 
and the updating rule for
(17) 
where
3.2 Complexity Analysis
For ease of analysis, we firstly define some symbols used in the following. Let and denote the total number of iterations and the lowest rank of the matrix , respectively. The size of is . The major computation cost of our proposed method contains two parts, calculating all ’s and updating . In terms of the formula (9) through (4) and (5), the computational complexity of Log algorithm is where is the number of terms in a discretized curves; therefore, the complexity of is at most and ’s computational complexity is . Thus the total for all the is . In each iteration of the Algorithm, the singular value thresholding is adopted to update the low rank matrix whose complexity is [12]. Suppose the algorithm is terminated after iterations, the overall computational complexity is given by
3.3 Convergence Analysis
Algorithm 1 is adopted from the algorithm proposed in [11]. However due to the terms of ’s in the objective function (11), the convergence theorem proved in [11] cannot be directly applied to this case as the linearization is implemented on both the augmented Lagrangian terms and the term involving ’s. Fortunately we can employ the revised approach, presented in [30], to prove the convergence for the algorithm. Without repeating all the details, we present the convergence theorem for Algorithm 1 as follows.
Theorem 1 (Convergence of Algorithm 1)
In all the experiments we have conducted, the algorithm converges very fast with .
4 Experiments
In this section we show three sets of experiments to evaluate the newly proposed cLRR. The performance of the proposed method is compared with the same type of subspace clustering algorithm LRR [12]. To compare segmentation accuracy we use the subspace clustering accuracy (SCA) metric [6], which is defined as
(18) 
Therefore a higher SCA means greater clustering accuracy.
The parameters used were fixed across all experiments with for LRR set at and for cLRR. A wide range of parameters were tested for each algorithm. Overall we found that the segmentation accuracy of LRR did not vary that much with changes in .
4.1 Synthetic Data
Mean  Median  Min  Max  

LRR  80.4%  83.33%  60%  91.67% 
CurveLRR  96.77%  98.33%  73.33%  100% 
To evaluate and confirm the effectiveness of the curve LRR method we first perform experimental evaluation using synthetic data. In this test three clusters were created consisting of twenty 1D curves of length 100. The curves in each cluster were sine waves, with each cluster corresponding to a unique frequency. Within each cluster progressive amounts of warping were applied. See Figure 1 for an example of data from three syntheticly generated clusters. Clustering was then performed on the data by applying curve LRR and segmenting the affinity matrix with nCUT. This experiment was repeated times with new data generated each time to obtain basic statistics. We compare against the baseline: LRR. Results are reported using subspace clustering accuracy and can be found in Table 1. Overall in this experiment Curve LRR outperforms conventional LRR by a significant margin.
4.2 Semisynthetic TIR Data
We assemble synthetic data from a library of pure infrared hyper spectral mineral data. For each cluster we pick one spectra sample from the library as a basis. Each curve basis is then randomly shifted and stretched in a random portion. This random warping is performed times to produce the curves for each cluster. See Figure 3 for an example of data used in this experiment. In this experiment we used three clusters. Again as in the previous experiment we repeated the test times. Results are reported in Table 2 and Figure 4.
The results show that LRR cannot accurately cluster data with this sort of nonlinear invariance, which is commonly found in this type of data due to impurities in the mineral samples. On the other hand cLRR perfectly clustered the data.
Mean  Median  Min  Max  

LRR  60.13%  60%  50%  71.67% 
CurveLRR  100%  100%  100%  100% 
4.3 Character Classification
In this experiment a collection of handwritten English characters were used to evaluate performance on a real world data set. The dataset consists of pen position data collected by a digitisation tablet at 200Hz, which is then converted to horizontal and vertical velocities [27, 26]. These 2D trajectory curves are normalised such that the mean of each curve is close to zero. See Figure 5 for some examples of this data. Figure 6 shows the example plots of curves used in the character classification experiment.
To evaluate performance twenty characters were randomly selected from three character classes. The data as originally released has been carefully produced and processed so that trajectories for each characters are extremely similar. Far more so than is realistic. For example the start time for each character has been aligned furthermore the writing speed, character size and variance in velocity over time is extremely consistent. Therefore to make the data more realistic we randomly globally shift each character so that their start times vary. Furthermore we randomly globally stretch and shrink each trajectory to account for different writing speeds, we also scale the trajectories by applying constant factors to account for character size and we lastly perform local warping (as done in the semisynthetic experiment) to account for variance in speed over time.
Since the data consists of multidimensional curves the X and Y trajectory curves were concatenated to form data usable for conventional LRR since it can only handle vectors. Results can be found in Table 3 and Figure 7. Once again, the cLRR clearly outperforms LRR in all metrics. Furthermore cLRR shows excellent performance with a median accuracy of over on an extremely challenging dataset.
Mean  Median  Min  Max  

LRR  52.33%  51.67%  43.33%  63.33% 
CurveLRR  86.33%  91.67%  70%  100% 
5 Conclusion
In this paper, we extended the conventional LRR model on Euclidean space to a new LRR model for the manifold of open curves. The new LRR formulation is based on the tangent space approximation to the manifold so that the classic data self expressive can be well preserved for the manifold of curves at relevant high accuracy. The resulting optimization problem can be solved using the LADMAP technique and algorithm convergence and complexity were presented. Finally we tested the new model by conducting experiments on synthetic, semisynthetic and real world data, and the experimental results show the outstanding performance against the conventional LRR. Our next work is further extended the LRR model to the manifold of general closed curves.
Acknowledgments
Funding information hidden for the review process.
References
 [1] P.A. Absil, R. Mahony, and R. Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, 2008.
 [2] M. T. Bahadori, D. Kale, Y. Fan, and Y. Liu. Functional subspace clustering with application to time series. In Proceedings of The 32nd International Conference on Machine Learning, pages 228–237, 2015.
 [3] C. Bishop. Pattern Recognition and Machine Learning. Information Science and Statistics. Springer, 2006.
 [4] J. F. Cai, E. J. Candès, and Z. Shen. A singular value thresholding algorithm for matrix completion. SIAM J. on Optimization, 20(4):1956–1982, 2008.
 [5] E. J. Candès, X. Li, Y. Ma, and J. Wright. Robust principal component analysis? Submitted for publication, Stanford University, 2010. http://wwwstat.stanford.edu/~candes/papers/RobustPCA.pdf.
 [6] E. Elhamifar and R. Vidal. Sparse subspace clustering: Algorithm, theory, and applications. IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(11):2765–2781, 2013.
 [7] F. Ferraty and Y. Romain, editors. The Oxford Handbook of Functional Data Analysis. Oxford University Press, 2011.
 [8] Y. Fu, J. Gao, X. Hong, and D. Tien. Low rank representation on Riemannian manifold of symmetrical positive definite matrices. In SIAM Conferences on Data Mining (SDM), pages 316–324, 2015.
 [9] S. H. Joshi, E. Klassen, A. Srivastava, and I. Jermyn. A novel representation for Riemannian analysis of elastic curves in . In IEEE Conference on Computer Vision and Pattern Recognition, pages 1–7, 2007.
 [10] Z. Lin, R. Liu, and H. Li. Linearized alternating direction method with parallel splitting and adaptive penalty for separable convex programs in machine learning. Machine Learning, 99:287–325, 2015.
 [11] Z. Lin, R. Liu, and Z. Su. Linearized alternating direction method with adaptive penalty for low rank representation. In Proceedings of NIPS, 2011.
 [12] G. Liu, Z. Lin, S. Yan, J. Sun, Y. Yu, and Y. Ma. Robust recovery of subspace structures by lowrank representation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(1):171–184, 2013.
 [13] H.G. Müller. International Encyclopedia of Statistical Science, chapter Functional data analysis, pages 554–555. Springer, 2011.
 [14] N. Parikh and S. Boyd. Proximal algorithms. Foundations and Trends in Optimization, 1(3):123–231, 2013.
 [15] F. Petitjean, G. Forestier, G. I. Webb, A. E. Nicholson, Y. Chen, and E. Keogh. Dynamic time warping averaging of time series allows faster and more accurate classification. in icdm, 2014. In International Conference on Data Mining, 2014.
 [16] T. Rakthanmanon. Addressing big data time series: Mining trillions of time series subsequences under dynamic time warping. ACM Transactions on Knowledge Discovery from Data, 7(3):1–31, 2013.
 [17] J. Ramsay and B. W. Silverman. Functional Data Analysis. Springer Series in Statistics. Springer, 2005.
 [18] D. Robinson. Functional analysis and partial matching in the square root velocity framework. PhD thesis, Florida State University, 2012.
 [19] J. Shi and J. Malik. Normalized cuts and image segmentation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 22:888–905, 2000.
 [20] A. Srivastava, E. Klassen, S. H. Joshi, and I. H. Jermyn. Shape analysis of elastic curves in Euclidean spaces. IEEE Transactionson Pattern Analysis and Machine Intelligence, 33(7):1415–1428, 2011.
 [21] A. Srivastava, W. Wu, S. Kurtek, E. Klassen, and J. S. Marron. Registration of functional data using FisherRao metric. varXiv:1103.3817, 2011.
 [22] J. Su and A. Srivastava. Rateinvariant analysus of trajectories on Riemannian manifolds with application in visual speech recognition. In Proceedings of International Conference on Computer Vision and Pattern Recognition, 2014.

[23]
J. Su, A. Srivastava, and F. W. Huffer.
Detection, classification and estimation of individual shapes in 2D and 3D point clouds.
Computational Statistics & Data Analysis, 58:227–241, 2013.  [24] J. D. Tucker, W. Wu, and A. Srivastava. Generative models for functional data using phase and amplitude separation. Computational Statistics and Data Analysis, 61:50–66, 2013.
 [25] B. Wang, Y. Hu, J. Gao, Y. Sun, and B. Yin. Low rank representation on grassmann manifolds: An extrinsic perspective. arXiv:1301.3529, 1:1–9, 2015.
 [26] B. Williams, M. Toussaint, and A. J. Storkey. Modelling motion primitives and their timing in biologically executed movements. In Advances in neural information processing systems, pages 1609–1616, 2008.
 [27] B. H. Williams, M. Toussaint, and A. J. Storkey. Extracting motion primitives from natural handwriting data. Springer, 2006.
 [28] L. Wu, A. Ganesh, B. Shi, Y. Matsushita, Y. Wang, and Y. Ma. Convex optimization based lowrank matrix completion and recovery for photometric stereo and factor classification. IEEE Transactions on Pattern Analysis and Machine Intelligence, XX:XXX–XXX, August 2012.
 [29] M. Yin, J. Gao, and Y. Guo. A nonlinear lowrank representation on Stiefel manifold. Electronics Letters, 51(10):749–751, 2015.
 [30] M. Yin, J. Gao, Z. Lin, Q. Shi, and Y. Guo. Graph dual regularized lowrank matrix approximation for data representation. IEEE Transactions on Image Processing, 24(12):4918–4933, 2015.
 [31] Z. Zhang, J. Su, E. Klassen, H. Le, and A. Srivastava. Videobased action recognition using rateinvariant analysis of covariance trajectories. arXiv:1503.06699v1, 1, 2015.