Year: 2020

Abstract:
We propose a novel local Subspace Indexing Model with Interpolation (SIM-I) for low-dimensional embedding of image data sets. Our SIM-I is constructed via two steps: in the first step we build a piece-wise linear affinity-aware subspace model under a given partition of the data set; in the second step we interpolate between several adjacent linear subspace models constructed previously using the “center of mass” calculation on Stiefel and Grassmann manifolds. The resulting subspace indexing model built by SIM-I is a globally non-linear low-dimensional embedding of the original data set. Furthermore, the interpolation step produces a “smoothed” version of the piece-wise linear embedding mapping constructed in the first step, and can be viewed as a regularization procedure. We provide experimental results validating the effectiveness of SIM-I, that improves PCA recovery for SIFT data set and nearest-neighbor classification success rates for MNIST and CIFAR-10 data sets.

heading: 1 INTRODUCTION
text: Subspace selection algorithms have been successful in many application problems related to dimension reduction (Zhou et al. (2010), Bian & Tao (2011), Si et al. (2010), Zhang et al. (2009)), with applications including, e.g., human face recognition (Fu & Huang (2008)), speech and gait recognition (Tao et al. (2007)), etc.. The classical approaches of subspace selection in dimension reduction include algorithms like Principle Component Analysis (PCA, see Jolliffe (2002)) and Linear Discriminant Analysis (LDA, see Belhumeur et al. (1997), Tao et al. (2009)). They are looking for globally linear subspace models. Therefore, they fail to estimate the nonlinearity of the intrinsic data manifold, and ignore the local variation of the data (Saul & Roweis (2003), Strassen (1969)). Consequently, these globally linear models are often ineffective for search problems on large scale image data sets. To resolve this difficulty, nonlinear algorithms such as kernel algorithms (Ham et al. (2004)) and manifold learning algorithms (Belkin et al. (2006), Guan et al. (2011)) are proposed. However, even though these nonlinear methods significantly improve the recognition performance, they face a serious computational challenge dealing with large-scale data sets due to the complexity of matrix decomposition at the size of the number of training samples.
Here we propose a simple method, Subspace Indexing Model with Interpolation (SIM-I), that produces from a given data set a piece-wise linear, locality-aware and globally nonlinear model of low-dimensional embedding. SIM-I is constructed via two steps: in the first step we build a piecewise linear affinity-aware subspace model under a given partition of the data set; in the second step we interpolate between several adjacent linear subspace models constructed previously using the “center of mass” calculation on Stiefel and Grassmann manifolds (Edleman et al. (1999), Kaneko et al. (2013), Marrinan et al. (2014)). The interpolation step outputs a “smoothed” version (Figure 1) of the original piece-wise linear model, and can be regarded as a regularization process. Compared to previously mentioned subspace methods, SIM-I enjoys the following advantages: (1) it captures the global nonlinearity and thus the local fluctuations of the data set; (2) it is computationally feasible to large-scale data sets since it avoids the complexity in matrix decomposition at the size of the number of training samples; (3) it includes a regularization step via interpolating between several adjacent pieces of subspace models. Numerical experiments on PCA recovery task for SIFT data set and classification tasks via nearest-neighbor method for MNIST and CIFAR-10 data sets further validate the effectiveness of SIM-I.

heading: 2 PIECE-WISE LINEAR LOCALITY PRESERVING PROJECTION (LPP) MODEL
text: If an image data point x ∈ RD is represented as a vector in a very high-dimensional space, then we want to find a low-dimensional embedding y = f(x) ∈ Rd, d << D such that the embedding function f retains some meaningful properties of the original image data set, ideally close to its intrinsic dimension. If we restrict ourselves to linear maps of the form y = WTx ∈ Rd, where the D × d projection matrix W = (wij)1≤i≤D,1≤j≤d (assuming full rank), then such a procedure is called a locally linear low-dimensional embedding (see Roweis & Saulm (2000); Van Der Maaten et al. (2009)). The target is to search for a “good” projection matrix W , such that the projection x 7→ y = WTx must preserve certain locality in the data set (this is called a Locality Preserving Projection, or LPP projection, see He & Niyogi (2003)). The locality is interpreted as a kind of intrinsic relative geometric relations between the data points in the original high-dimensional space, usually represented by the affinity matrix S = (sij)1≤i,j≤n (which is a symmetric matrix with non-negative terms). As an example, given unlabelled data points x1, ..., xn ∈ RD, we can take sij = exp ( −‖xi−xj‖ 2
2σ2
) when ‖xi − xj‖ < ε and sij = 0 otherwise. Here σ > 0 and ε > 0 is a
small threshold parameter, and ‖xi−xj‖ is the Euclidean norm in RD. Based on the affinity matrix S = (sij), the search for the projection matrix W can be formulated as the following optimization problem
min W
φ(W ) = 1
2 n∑ i,j=1 sij‖yi − yj‖2 , (1)
in which yi = WTxi and yj = WTxj and the norm ‖yi − yj‖ is taken in the projected space Rd. Usually when ‖xi − xj‖ is large, the affinity sij will be small, and vice versa. Thus (1) is seeking for the embedding matrix W such that close pairs of image points xi and xj will be mapped to close pairs of embeddings yi = WTxi and yj = WTxj , and vice versa. This helps to preserve the local geometry of the data set, a.k.a the locality. To solve (1), we introduce a weighted fully-connected graph G where the vertex set consists of all data points x1, ..., xn and the weight on the edge connecting xi and xj is given by sij ≥ 0. Consider the diagonal matrix D = diag(D11, ...,Dnn) where
Dii = n∑ j=1 sij , and we then introduce the graph Laplacian L = D − S. Then the minimization
problem (1), together with the normalization constraint n∑ i=1 Diiy2i = 1, reduces to the following generalized eigenvalue problem (see He & Niyogi (2003))
XLXTw = λXDXTw , (2) where X = [x1, ..., xn] ∈ RD×n. Assume we have obtained an increasing family of eigenvalues 0 = λ0 < λ1 ≤ ... ≤ λn−1. Let the corresponding eigenvectors be w0, w1, ..., wn−1. Then the low-dimensional embedding matrix can
be taken as W = [w1, ..., wd] (see He et al. (2005)). By choosing different affinity matrices S = (sij), the above LPP framework includes many commonly seen practical examples. For example, if the data x1, ..., xn are not labelled, then we can take sij = 1
n and (2) produces the classical Principle
Component Analysis (PCA). For labelled data forming subsets X1, ...,Xm with same labels in each subset, we can take sij = 1
nk when xi, xj ∈ Xk, and sij = 0 other-wise. Here nk is the cardinality
of Xk. This will produce Linear Discriminant Analysis (LDA). The detailed justifications of these connections can be found in He et al. (2005).
Given an input data set X = {x1, ..., xn} where each xi ∈ RD, either labelled or unlabelled, we can apply a k-d tree (Bentley (1975), Wang et al. (2011)) based partition scheme to divide the whole data set X into non-overlapping subsets C1, ..., C2h where h is the depth of the tree. Conventional subspace selection algorithms could be applied on the whole sample space before the whole space is partitioned and indexed. For example, we can first apply a PCA to X , which selects the first d bases [a1, ..., ad] with largest variance. Based on these bases, the covariance information obtained from global PCA is utilized in the indexing as follows: (1) we project all sample points x1, ..., xn onto the maximum variance basis a1, then we find the median value m1 of the projected samples, and split the whole collection of data along a1 at m1, i.e., split the current node into left and right children; (2) starting from level i = 2, for each left and right child, project the whole collection of data along the i-th maximum variance basis ai, find the median value mi, and split all the children at mi; (3) increase the level from i to i+ 1 and repeat (2) until i = h reaches the bottom of the tree. We collect all the 2h children at level i = h and obtain the disjoint subsets C1, ..., C2h . Each subset Ck, k = 1, 2, ..., 2h consists of a family of input data in RD. Based on them, using the above LPP framework, for each Ck, a low-dimensional embedding matrix Wk ∈ RD×d can be constructed. In this way, over the whole data set X , we have constructed a piece-wise linear low-dimensional embedding model f(x) : RD → Rd (see Figure 1(a)) where x ∈ X . This model is given by the linear embedding matrices W1, ...,W2h ∈ RD×d. The above model construction can be regarded as a training process from the data set X . For a given test data point x ∈ RD, not included in X , we can find the closest subset Ck(x) to it, by selecting the index k = k(x) ∈ {1, ..., 2h} with the smallest distance ‖x−mk(x)‖. Here mk is the mean among all data points in the subset Ck. With the subset Ck(x) chosen, we map the test point x ∈ RD to its low-dimensional embedding f(x) = WTk(x)x ∈ R
d. Such a procedure extended the piece-wise linear embedding model f(x) : RD → Rd to all testing data points in RD.

heading: 3 CALCULATING THE “CENTER OF MASS” ON STIEFEL AND GRASSMANN MANIFOLDS
text: Subspace Indexing (Wang et al. (2011)) provides a d-dimensional representation of the data set {x1, ..., xn} by the subspace span(w1, ..., wd) = {WTx, x ∈ RD} generated from the linear embedding matrix W ∈ RD×d. In this case, we are only interested in the column space of W , so we can assume that w0, w1, ..., wn−1 is an orthonormal basis 1. Such a matrix W belongs to the Stiefel manifold, defined by
Definition 1 (Stiefel manifold) The compact Stiefel manifold St(d,D) is a submanifold of the Eucilidean space RD×d such that
St(d,D) = {X ∈ RD×d : XTX = Id} . (3)
As an example, if we are interested in signal recovery using low-dimensional PCA embedding, the projections we calculated from PCA analysis will be on Stiefel manifolds. However, for classification tasks, the exact distance information is less important than label information. In this case, two such Stiefel matricesW1 andW2 produce the same embedding ifW1 = W2Od for someOd ∈ O(d), whereO(d) is the group of orthogonal matrices in dimension d. In this case, the relevant embedding we obtained is a point on the Grassmann manifold, defined by
1If this is not the case, we can replace the matrix [w0 w1 ... wn−1] by the Q matrix of the QR-decomposition of itself, without changing the corresponding subspace.
Definition 2 (Grassmann manifold) The Grassmann manifold Gr(d,D) is defined to be the quotient manifold Gr(d,D) = St(d,D)/O(d). A point on Gr(d,D) is defined by an equivalence class [W ] = {WOd, Od ∈ O(d)} where W ∈ St(d,D).
Given a family of elements on the Stiefel or Grassmann manifold, the center-of-mass is defined as an element on the same manifold that minimizes the functional given by the weighted sum of square distances. To be precise, we have
Definition 3 (Stiefel and Grassmann center-of-masses) Given a sequence of matrices W1, ...,Wl ∈ St(d,D) and a sequence of weights w1, ..., wl > 0, the Stiefel center-of-mass with respect to the distance d(W1,W2) on St(d,D) is defined as a matrix Wc = W Stc (W1, ...,Wl;w1, ..., wl) ∈ St(d,D) such that
Wc = W St c (W1, ...,Wl;w1, ..., wl) ≡ arg min
W∈St(d,D) l∑ j=1 wjd 2(W,Wj) . (4)
Similarly, if the corresponding equivalent classes are [W1], ..., [Wl] ∈ Gr(d,D), then the Grassmann center-of-mass with respect to the distance d([W1], [W2]) on Gr(d,D) is defined as the equivalence class [Wc], where Wc = WGrc (W1, ...,Wl;w1, ..., wl) ∈ St(d,D) is such that
Wc = W Gr c (W1, ...,Wl;w1, ..., wl) ≡ arg min
W∈St(d,D) l∑ j=1 wjd 2([W ], [Wj ]) . (5)
The distances d(W1,W2) or d([W1], [W2]) can be taken in different ways. For example, for W1,W2 ∈ St(d,D), one way is to consider d(W1,W2) = dF (W1,W2) = ‖W1 −W2‖F , the matrix Frobenius norm of W1 −W2. One can also take a more intrinsic distance, such as the geodesic distance between W1 and W2 on the manifold St(d,D) with the metric given by embedded geometry (see Edleman et al. (1999)). For [W1], [W2] ∈ Gr(d,D), one way is to consider the projected Frobenius norm d([W1], [W2]) = dpF ([W1], [W2]) = 2−1/2‖W1WT1 −W2WT2 ‖F . There are also many other choices, such as using the principle angles between the subspaces, chordal norms, or other types of Frobenius norms (see Edleman et al. (1999, Section 4.3)).
With respect to matrix Frobenius norm and projected Frobenius norm, the Stiefel and Grassmann center-of-masses can be calculated explicitly in the following theorems.
Theorem 1 (Stiefel center-of-mass with respect to Frobenius norm) We consider the singular value decomposition of the matrix ∑l j=1 wjWj = O1∆O2, where O1 ∈ O(D) and O2 ∈ O(d),
∆ =
( diag(λ1, ..., λd)d×d
0(D−d)×d
) and λ1 ≥ ... ≥ λd ≥ 0 are the singular values. Then the Stiefel
center-of-mass with respect to the distance given by Frobenius norm d(W1,W2) = ‖W1−W2‖F is given by Wc = O1ΛO2 where Λ = (
diag(1, ..., 1)d×d 0(D−d)×d
) .
Theorem 2 (Grassmann center-of-mass with respect to projected Frobenius norm) Set Ωj =(∑l j=1 wj )−1 wj . We consider the singular value decomposition of the symmetric matrix
l∑ j=1 ΩjWjW T j = Q∆Q T where Q ∈ O(D) and ∆ = diag(σ21 , ..., σ2D), σ21 ≥ ... ≥ σ2D ≥ 0. Then the Grassmann center-of-mass with respect to the distance given by projected Frobenius norm dpF ([W1], [W2]) = 2−1/2‖W1WT1 − W2WT2 ‖F is the equivalence class [Wc] determined
by Wc = QΛ, where Λ = (
diag(1, ..., 1)d×d 0(D−d)×d
) .
4 SIM-I: INTERPOLATING THE LPP MODEL FAMILY
Recall that we have developed a piece-wise linear embedding model f(x) : RD → Rd over the data set X = {x1, ..., xn}. The embedding f(x) corresponds to a family of subspace indexing models
W1, ...,W2h ∈ St(d,D) (e.g. for PCA signal recovery tasks) or [W1], ..., [W2h ] ∈ Gr(d,D) (e.g. for classification tasks). Each subspace model Wk is built from LPP embedding using the subset Ck ⊂ X developed from k-d tree and h is the depth tree. Given a test point x ∈ RD that does not lie in X , we can map it to the low-dimensional embedding f(x) = WTk(x)x ∈ R
d. The index k(x) corresponds to the subset Ck(x) that lie closest to x. In practice, we can first compute the means mk over all the data points in the subset Ck for each k = 1, 2, ..., 2h and sort the distances ‖x−mk‖ in ascending order ‖x−mk1(x)‖ ≤ ... ≤ ‖x−mk2h (x)‖, {k1(x), ..., k2h(x)} = {1, ..., 2
h}. We then take k(x) = k1(x) to be the index k corresponding to the shortest distance. This is effective when the test point x lies significantly close to one of the subsets Ck(x), see Figure 2(a).
Algorithm 1 SIM-I: Subspace Indexing Model with Interpolation 1: Input: Data set X = {x1, ..., xn ∈ RD} and its corresponding affinity matrix S =
(sij)1≤i,j≤n; test point x ∈ RD; threshold ratio rthr > 1; tree depth h; parameter K > 0
2: Using an initial PCA and a k-d tree based partition scheme, decompose the data set X into subsets C1, ..., C2h , where h is the depth of the tree 3: For each subset Ck, calculate its mean (center) mk ∈ RD and its LPP embedding matrix Wk ∈ St(d,D) based on the affinity matrix S 4: Sort the distances ‖x − mk‖ in ascending order ‖x − mk1(x)‖ ≤ ... ≤ ‖x − mk2h (x)‖, {k1(x), ..., k2h(x)} = {1, ..., 2h} 5: Determine I , which is the first sub-index i of ki(x) such that ‖x−mkI+1(x)‖ > rthr‖x−mk1(x)‖
6: Set ji(x) = ki(x) for i = 1, 2, ..., I and obtain the embedding matrices Wj1(x), ...,WjI(x) ∈ St(d,D) or their corresponding subspaces [Wj1(x)], ..., [WjI(x)] ∈ Gr(d,D), together with the weights wi = exp(−K‖x−mji(x)‖2) > 0 for i = 1, ..., I 7: Find a center-of-mass Wc = W Stc (Wj1(x), ...,WjI(x);w1, ..., wI) (Stiefel case) or [Wc] = [WGrc (Wj1(x), ...,WjI(x);w1, ..., wI)] (Grassmann case) according to Definition 3 and Theorems 1 and 2. 8: Output: The low-dimensional embedding f(x) = WTc x ∈ Rd
However, for a general test point x ∈ RD, it might happen that this point lies at approximately the same distances to the centers of each of the several different subsets adjacent to x (see Figure 2(b)).
In this case, we aim to interpolate between several subspace indexing modelsWj1(x), ...,WjI(x). To do this, we first find the subspace indexes j1(x), ..., jI(x) from the first I subsets Cj1(x), ..., CjI(x) closest to x, i.e., j1(x) = k1(x), ..., jI(x) = kI(x) given the sorted distances ‖x−mk‖ mentioned above. In practice, the number I = I(x) is depending on x and can be chosen in the following way: I is the first sub-index i of ki(x) such that ‖x−mkI+1(x)‖ > rthr‖x−mk1(x)‖, where rthr > 1 is a threshold ratio that can be tuned. We then pick the weights aswi = exp(−K‖x−mji(x)‖2) for some K > 0 and i = 1, 2, ..., I . This is indicating that the closer x is to Cji(x), the heavier weights we assign to Wji(x) in the interpolation process. Given the embedding matrices Wj1(x), ...,WjI(x) ∈ St(d,D) or their corresponding subspaces [Wj1(x)], ..., [WjI(x)] ∈ Gr(d,D), together with the weights w1, ..., wI > 0, we find a center-of-mass Wc = W Stc (Wj1(x), ...,WjI(x);w1, ..., wI) (Stiefel case) or [Wc] = [WGrc (Wj1(x), ...,WjI(x);w1, ..., wI)] (Grassmann case) according to Definition 3 and Theorems 1 and 2. Finally, we map the test point x to the low-dimensional embedding f(x) = WTc x ∈ Rd. Notice that when I = 1, the interpolation procedure reduces to projecting x using Wk1(x) calculated from LPP analysis on the closest subset only. In general, the whole interpolation procedure can be regarded as providing a regularized version of the piece-wise linear embedding we discussed in Section 2 (also see Figure 1). We summarize our interpolation method as the SIM-I Algorithm 1.

heading: 5 EXPERIMENTS

heading: 5.1 PCA RECOVERY FOR SIFT DATA SET
text: SIFT (Scale Invariant Feature Transform, see Lowe (2004)) data set is a data set that computes for each keypoint a real-valued descriptor, based on the content of the surrounding patch in terms of local intensity gradients. Given its remarkable performance, SIFT has been often used as starting point for the creation of other descriptors. The final SIFT descriptor consists of 128 elements. This means that each data point in SIFT data set has dimension D = 128, and we pick the embedding dimension d = 16. The original SIFT data set has 10068850 data samples, that form the data set sift sample. We randomly collect ntrain = 200 × 213 elements from these sample points as our target data set X = sift train = {x1, ..., xntrain ∈ RD}. We consider the recovery efficiency of PCA embedding of the SIFT data set. Let x be a point in sift sample and let W be a Stiefel matrix in St(16, 128). The projection of x onto the 16-dimensional space is then denoted by y = WTx. By recovery we meant to consider the point x̂ = (WT )−y where (WT )− is the Moore-Penrose pseudo-inverse of WT . The recovery efficiency can then be quantified by the recovery error ‖x− x̂‖, where the Euclidean norm is computed in R128. We pick h = 13 so that 2h = 8192. Then we decompose sift train into 8192 subsets C1, ..., C8192 using k-d tree based partition. For each subset Ck, we calculate the mean mk ∈ R128 and we obtain a PCA embedding matrix Wk ∈ St(16, 128). We sort the distances ‖x −mk‖, k = 1, 2, ..., 8192 in ascending order so that ‖x −mk1‖ ≤ ... ≤ ‖x −mk8192‖, where {k1, ..., k8192} = {1, ..., 8192}. Then we find among the subset means mk, k = 1, 2, ..., 8192 the first I nearest to x, with their indexes denoted by ki = ki(x), i = 1, .., I . The number I = I(x) is depending on x and is chosen in the following way: I is the first sub-index i of ki such that ‖x−mkI+1‖ > rthr‖x−mk1‖, where rthr > 1 is a threshold ratio. We pick rthr = 2. For test data set, we randomly pick from sift sample\sift train a subset of size ntest = 500, and we denote the data set as sift test. For each test point x ∈ sift test, we consider sending it to the nearest subset Ck1 and the corresponding Stiefel matrix is Wk1 ∈ St(16, 128). We can then consider the recovery point x̂ = (WTk1)
−WTk1x and the benchmark recovery error Error bm = ‖x− x̂‖. Consider the alternative recovery scheme using our method SIM-I. We calculate the weights wj = exp(−K‖x − mkj‖2) for j = 1, ..., I and we choose the constant K = 10−8. We then find the Stiefel center-of-mass Wc = W Stc (Wk1(x), ...,WkI(x);w1, ..., wI) using Theorem 1 and taking the distance function to be the matrix Frobenius norm d(W1,W2) = ‖W1 −W2‖F . We consider the recovery point x̂c = ((Wc)T )−(Wc)Tx and the recovery error Error c = ‖x− x̂c‖. Over the test set sift test, we find that for about 94.2% test points, Error c < Error bm, which implies that SIM-I improved the efficiency of recovery. The empirical average Error c is
402.089506, and the empirical average Error bm is 454.452314. Figure 3 plots the Error c and Error bm (vertical axis) as functions of the test sample indexes. The red curve is for Error c and blue curve is for Error bm, where we have sorted Error bm in ascending order and reordered the test indexes correspondingly. Figure 4 gives the differences Error c − Error bm in descending order. It can be apparently seen that Error c− Error bm < 0 for most of test samples.

heading: 5.2 NEAREST-NEIGHBOR CLASSIFICATION FOR MNIST AND CIFAR-10
text: Here we consider a labeled data set data = {(xi, yi), i = 1, 2, ..., N} where xi ∈ RD are the inputs and yi ∈ N are the labels. We randomly select the training set data train = {(xi, Yi), i = 1, 2, ..., ntrain} and the test set data test = {(ai, bi), i = 1, 2, ..., ntest}, and we make them disjoint. We first project the training set onto a kd PCA-dimensional subspace via a standard PCA.
Based on this initial embedding, using a kd-tree with height h, we divide data train into 2h clusters Ck, k = 1, 2, ..., 2h. For each Ck, we find an LPP embedding matrixWk ∈ St(kd LPP, D) by setting the affinity matrix to be sij = exp(−‖xi− xj‖2) if xi, xj ∈ Ck are in the same class and sij = 0 otherwise. Since we are having a classification problem, we can identify each Wk by the subspace it spans, i.e., we consider [Wk] which is the equivalence class of Wk in Gr(kd LPP, D).
As before, we compute in RD the means mk of each cluster Ck. For a test point x ∈ data test, we sort the distances ‖x − mk‖, k = 1, 2, ..., 2h in ascending order so that ‖x − mk1‖ ≤ ... ≤ ‖x −mk
2h ‖, where {k1, ..., k2h} = {1, ..., 2h}. Then we find among the cluster means mk, k =
1, 2, ..., 2h the first I nearest to x, with their indexes denoted by ki = ki(x), i = 1, ..., I . The number I = I(x) is depending on x and is chosen in the following way: I is the first sub-index i of ki such that ‖x −mkI+1‖ > rthr‖x −mk1‖, where rthr > 1 is a threshold ratio. The parameter rthr can be treated as a hyper-parameter that we can tune here.
For baseline method, we do a nearest-neighbor classification for x on a low-dimensional embedding of the cluster Ck1 , and we pick the number of nearest-neighbors to be knn ≥ 1. Indeed we project x and all training data points in Ck1 using Wk1 , and perform nearest-neighbor classification on the resulting projection.
For our method SIM-I, we take the union Ck1 ∪ ... ∪ CkI . Recall that each Cki corresponds to an LPP embedding projection matrix Wki . We set the weights wi = exp(−K‖x − mki‖2), and we pick K = 10−8. We compute a center-of-mass Wc of the projection matrices Wki with weights wi, i = 1, 2, ..., I using the Grassmann center-of-mass method, where the distance is taken as the projected Frobenius norm of Grassmann matrices, i.e., d(W1,W2) = 2−1/2‖W1WT1 −W2WT2 ‖F and ‖ • ‖F is the matrix Frobenius norm. We obtain Wc = WGrc (Wk1(x), ...,WkI(x);w1, ..., wI) and we project x and all training data points in the union Ck1 ∪ ... ∪ CkI using Wc. We then perform a nearest-neighbor classification for x on this low-dimensional embedding with the number of nearest-neighbors being equal to knn.
We have been experimenting on 2 different data sets, for both of them we pick kd PCA = 128, kd LPP = 100: (1) The MNIST data set, with N = 70000, ntrain = 60000, ntest = 10000 data points and D = 784; (2) The CIFAR-10 data set, with N = 60000, ntrain = 50000, ntest = 10000 data points and D = 3072. Table 1 shows the results, where the last 2 columns are the nearestneighbor classification success rates for baseline method, and for SIM-I using Grassmann centerof-mass and the rows are for different experiments. The first 4 columns are for data sets, kd-tree height, threshold value rthr, the number of nearest-neighbors knn, respectively. We slightly tuned rthr to reach best performances. Clearly, SIM-I has its advantage.

heading: A PROOF OF THEOREM 1
text: Let D ≥ d ≥ 1. Recall St(d,D) stands for the Stiefel manifold, that is, each matrix in St(d,D) is a D by d matrix with columns being orthogonally normal. For any matrix M , recall that ‖M‖2F = [tr(MTM)]1/2 is the Frobenius norm of M .
Theorem 1 (Stiefel center-of-mass with respect to Frobenius norm) We consider the singular value decomposition of the matrix l∑
j=1
wjWj = O1∆O2, where O1 ∈ O(D) and O2 ∈ O(d),
∆ =
( diag(λ1, ..., λd)d×d
0(D−d)×d
) and λ1 ≥ ... ≥ λd ≥ 0 are the singular values. Then the Stiefel
center-of-mass with respect to the distance given by Frobenius norm d(W1,W2) = ‖W1−W2‖F is given by Wc = O1ΛO2 where Λ = (
diag(1, ..., 1)d×d 0(D−d)×d
) .
Proof. For W ∈ St(d,D) we define f(W ) = l∑
j=1
wj‖W − Wj‖2F , and we are looking for a
minimizer Wc of f on St(d,D). Write
‖W −Wj‖2F = tr[(W −Wj)T (W −Wj)] = tr(WTW ) + tr(WTj Wj)− 2tr(WTWj) = 2p− 2 tr(WTWj).
It follows that
f(W ) = 2p l∑
j=1
wj − 2 l∑
j=1
wj tr(WTWj)
= 2p l∑
j=1
wj − 2 tr(WTB)
where B := l∑
j=1
wjWj . By singular value decomposition, there are D × D orthogonal matrix O1 and d× d orthogonal matrix O2 such that B = O1∆O2, where ∆ := (
diag(λ1, ..., λd)d×d 0(D−d)×d
) , and
λ1 ≥ · · · ≥ λd ≥ 0 are the singular values of B. It follows that tr(WTB) = tr(∆O2WTO1). Observe that (O2WTO1)T ∈ St(d,D). Write O2WTO1 = (cij)d×D. Then D∑ j=1 c2ij = 1 for each i. Hence
tr(WTB) = d∑ i=1 λicii ≤ d∑ i=1 λi ,
with equality holds when O2WTO1 = (Id, O) where O is the d × (D − d) matrix with all entries being zero. This says that Wc = O1ΛO2 where Λ = (
diag(1, ..., 1)d×d 0(D−d)×d
) is the maximizer of
tr(WTB). The so obtained Wc serves as the minimizer of f(W ) on St(d,D) and is thus the centerof-mass.

heading: B PROOF OF THEOREM 2
text: Recall that Gr(d,D) = St(d,D)/O(d). Every point W on St(d,D) will correspond to an equivalence class [W ] = {WOd : Od ∈ O(d)} which is a point on Gr(d,D). To represent points on Gr(d,D) using matrices, we notice that every point on Gr(d,D) corresponds to a unique choice of matrix WWT where W ∈ St(d,D). In this way, we can define the projected Frobenius distance between two classes [W1] and [W2] in Gr(d,D) as d2pF ([W1], [W2]) = 2 −1/2‖W1WT1 −W2WT2 ‖F .
Theorem 2 (Grassmann center-of-mass with respect to projected Frobenius norm) Set Ωj =(∑l j=1 wj )−1 wj . We consider the singular value decomposition of the symmetric matrix
l∑ j=1 ΩjWjW T j = Q∆Q T where Q ∈ O(D) and ∆ = diag(σ21 , ..., σ2D), σ21 ≥ ... ≥ σ2D ≥ 0. Then the Grassmann center-of-mass with respect to the distance given by projected Frobenius norm dpF ([W1], [W2]) = 2−1/2‖W1WT1 − W2WT2 ‖F is the equivalence class [Wc] determined
by Wc = QΛ, where Λ = (
diag(1, ..., 1)d×d 0(D−d)×d
) .
Proof. Set M = WWT and Mj = WjWTj , then we are looking for a minimizer Wc ∈ St(d,D) of f(W ) ≡ ∑l j=1 wjd 2 pF ([W ], [Wj ]) = 2 −1 l∑
j=1
wj‖M −Mj‖2F . For two D ×D matrices M1,M2
we define 〈M1,M2〉 = tr(MT1 M2). It is easy to verify that 〈, 〉 is an inner product and ‖M‖2F = 〈M,M〉. Thus we write
2f(W ) = l∑
j=1
wj‖M −Mj‖2F
= l∑
j=1
wj〈M,M〉 − 2 l∑
j=1
wj〈M,Mj〉+ l∑
j=1
wj〈Mj ,Mj〉
= l∑
j=1
wj
[ 〈M,M〉 − 〈M, 2
l∑ j=1 ΩjMj〉
] +
l∑ j=1 wj〈Mj ,Mj〉
= l∑
j=1
wj〈M − l∑
j=1
ΩjMj ,M − l∑
j=1
ΩjMj〉 − l∑
j=1
wj ∥∥∥∥∥ l∑j=1 ΩjMj ∥∥∥∥∥ 2
F
+ l∑
j=1
wj‖Mj‖2F
= l∑
j=1
wj ∥∥∥∥∥M − l∑j=1 ΩjMj ∥∥∥∥∥ 2
F
− m∑ j=1 wj ∥∥∥∥∥ l∑j=1 ΩjMj ∥∥∥∥∥ 2
F
+ l∑
j=1
wj‖Mj‖2F .
So the minimum is taken when ∥∥∥∥∥M − l∑j=1 ΩjMj ∥∥∥∥∥ 2
F
is minimized. Since eachWjWTj is a symmetric
matrix, we consider the SVD decomposition l∑
j=1
ΩjWjW T j = Q∆Q T where Q ∈ O(D) and ∆ =
diag(σ21 , ..., σ 2 D), σ 2 1 ≥ σ22 ≥ ... ≥ σ2D ≥ 0. Thus it suffices to set Wc = arg min W∈St(d,D) ‖WWT − Q∆QT ‖2F . Since W ∈ St(d,D), we have WWT = PV PT where P is an orthogonal matrix of size D ×D and V = diag(1, 1, .., 1, 0, 0, ..., 0) is an D ×D matrix with rank(V ) = d. Moreover, P can be chosen as P = (W Z) where Z is an D × (D − d) matrix. So
min W∈St(d,D) ‖WWT −QV QT ‖2F = min P∈O(D) ‖PV PT −Q∆QT ‖2F .
Let the orthgonal matrix O = Q−1P . Then we further have
min P∈O(n) ‖PV PT −Q∆QT ‖2F = min O∈O(n) ‖Q(OV OT −∆)QT ‖2F = min O∈O(n) tr(Q(OV OT −∆)2QT )
= min O∈O(n)
tr(OV OT −∆)2 .
We then show that min O∈O(n) tr(OV OT −∆)2 is achieved at the orthogonal matrix
O∗ =
( Op×p 0
0 O(n−p)×(n−p)
)
for any Op×p ∈ O(p) and any O(n−p)×(n−p) ∈ O(n − p). Indeed, we have tr(OV OT − ∆)2 = tr(OV OT ) + tr(∆2)− 2tr(∆OV OT ), so that we only have to maximize tr(∆OV OT ) = D∑ i=1 σ2i cii
where OV OT = (cij)D×D. Let O = (vij)D×D, then it is easy to calculate that cii = d∑ j=1 v2ij .
Moreover since D∑ j=1 v2ij = 1 for all i and D∑ i=1 v2ij = 1 for all j, we know that 0 ≤ cii ≤ 1 and
D∑ i=1 cii = d. Thus from σ21 ≥ ... ≥ σ2D ≥ 0 we see tr(∆OV OT ) = D∑ i=1 σ2i cii ≤ d∑ i=1 σ2i with equality if and only if cii = 1 for i = 1, ..., d and cii = 0 for i = d + 1, ..., D. This gives that, for
O = (vij)D×D, we have d∑ j=1 v2ij = 1 when i = 1, ..., d and D∑ j=d+1 v2ij = 1 when i = d + 1, ..., D.
Thus the minimum of min O∈O(n)
tr(OV OT −∆)2 is achieved at O∗ = ( Od×d 0
0 O(D−d)×(D−d)
) for
any Od×d ∈ O(d) and any O(D−d)×(D−d) ∈ O(D − d). We can pick Od×d = diag(1, ..., 1)d×d, and since P ∗ = QO∗ = (Wc Z) where Z is an D × (D − d) matrix, we get the conclusion of the Theorem about Wc.

Review summaries:
This paper proposes a subspace indexing model with interpolation (SIM-I) for dimension reduction. To capture the global nonlinearity and the local variation of data, SIM-I split the global space into a collection of local disjoint partitions by a kd-tree style scheme and build a subspace indexing model (e.g., LPP) for each local partition. The main contribution of this paper is to propose an interpolation technique to combine several subspace indexing models to determine a more promising model for low-dimensional embedding. Two experiments on PCA recovery and nearest neighbour classification validate the effectiveness of SIM-I.

This paper presents an approach to low-dimensional embedding of data, with an emphasis on image datasets. The algorithm proceeds by first finding several local subspaces that fit the data and then tying these subspaces together using the center of mass calculation on the Stiefel/Grassmann manifold. The method appears to be computationally efficient, since it can be applied to fairly large image datasets, though no formal analysis is given. The main drawback of the paper is that the empirical results compare only to PCA. Several more modern methods for low-dimensional embeddings exist (e.g., T-SNE, UMAP), and it would be appropriate to compare to these. UMAP also allows for the incorporation of labeled data.

The authors propose an approach towards computing  globally non-linear, low-dimensional embeddings of high dimensional data. In particular, they consider a sub-space indexing model with interpolation (called SIM-I) which consists of two steps. First, a locally linear model of the data is generated (using kD-trees and PCA) and second, averaging operations on the respective projection matrices are carried out. These are computed based on center of mass computations on Stiefel or Grasammanian manifolds. Anectdotal experimental evidence shows that this approach outperforms simple baseline techniques.