Fast Classical Scaling
Abstract
Multidimensionalscaling (MDS) is a dimensionality reduction tool used for information analysis, data visualization and manifold learning. Most MDS procedures find embedding of data points in low dimensional Euclidean (flat) domains, such that distances between the points are as close as possible to given interpoints dissimilarities. We present an efficient solver for Classical Scaling, a specific MDS model, by extrapolating the information provided by distances measured from a subset of the points to the rest. The computational and space complexities of the new MDS procedures are be thereby reduced from quadratic to quasilinear in the number of data points. Incorporating both local and global information about the data allows us to construct a low rank approximation to the intergeodesic distances between the data points. As a byproduct, the proposed method allows for efficient computation of geodesic distances. Finally, we show how to apply our method to two geometric analysis applications and obtain state of the art results.
1 Introduction
need for data simplification and dimensionality reduction tools constantly expands with the increasing availability of digital information. SelfOrganizing Map (SOM) [17], and Local Coordinate Coding [34, 35] are examples of such data reduction techniques, that can be used to simplify the data and hence reduce computational and space complexities of procedures dealing with it. Multidimensional Scaling (MDS) [5] is an example of such methods that embeds the data into a low dimensional Euclidean space, while attempting to preserve the distance between each pair of points. Flat embedding is a fundamental step in various applications in the fields of data mining, statistics, manifold learning, nonrigid shape analysis, and more. For example, in [21] Panazzo et al. demonstrated how to generalize weighted averages to surfaces using MDS. It enabled them to efficiently compute splines on surfaces, solve remeshing problems, find shape correspondences, transfer texture and evaluate morphing. Zigelman et al. [36] used MDS for texture mapping, while Schwartz et al. [27, 11], utilized it to flatten models of monkeys’ cortical surfaces for further analysis. In [28, 25, 23], MDS was applied to image and video analysis. In [12], MDS was used to embed surfaces into a Euclidean space, such that the Euclidean distance between each pair of embedded points is as similar as possible to the geodesic distance between the corresponding surface points. This embedding is referred as a canonical form and is invariant to isometric deformations of the surface. Such canonical forms reveal the intrinsic structure of the surface, and can be used to simplify tasks such as nonrigid objects matching and classification.
The bottleneck step in MDS is computing and storing all pairwise distances. For example, when dealing with geodesic distances computed on surfaces, this step is time consuming and in some cases impractical. The fast marching method [16] is an example of one efficient procedure for approximating the geodesic distance map between all pairs of points on a surface. Its complexity is when used in a straightforward manner. Even when using such an efficient procedure, the complexities involved in computing and storing these distances are at least quadratic in the number of points.
One way to reduce these complexities is by considering only local distances between nearby data points. Such attempts were made, for example, in Locally Linear Embedding (LLE, [24]) Hessian Locally Linear Embedding (HLLE, [9]), and in Laplacian eigenmaps [3]. This way, the effective time and space complexities can be reduced to . However, only the local structure is captured and the global geometry is ignored. In [31], De Silva and Tenenbaum suggested to compute only distances between a small subset of landmark points. The landmarks were then embedded ignoring of the rest of the distances, and then the rest of the points were interpolated into the target flat space.
The Nyström method [2] is an efficient technique to reconstruct a positive semidefinite matrix using only a small subset of its columns and rows. In [33, 7] only a few columns of the pairwise distance matrix were computed and then the Nyström method was used to interpolate the rest of the matrix. The main differences between these methods are the columns sampling schemes and the way to compute the initial geodesic distances. In addition, a regularization term was used in [7] for the pseudoinverse in the Nyström method, which further improves the reconstruction. [19] combined the Nyström method with KernelPCA [26] and demonstrated an efficient application of mesh segmentation. Recently, Spectral MDS (SMDS) [1] obtained stateoftheart results for efficiently approximating the embedding of MDS. In this method, complexities are reduced by translating the problem into the spectral domain, and representing all intergeodesic distances considering only the first eigenvectors of the LaplaceBeltrami operator.
We propose two methods for improving both accuracy and time complexities of MDS when compared to existing methods. Inspired by some ideas of SMDS, FastMDS (FMDS) [29] interpolates the distance map from a small subset of landmarks using a smoothness assumption as a prior. As opposed to the SMDS, the problem is solved in the spatial domain and accuracy is improved since no eigenvectors are omitted. NyströmMDS (NMDS) [30] further improves the solution by learning the distance interpolation coefficients from an initial set of computed distance. In this paper, similar to most above mentioned MDS acceleration techniques, Classical Scaling is used as our MDS model. It provides an optimal closed form solution. Both methods explored in this paper reconstruct the full distance matrix as a lowrank multiplication of small matrices, and reformulate classical scaling through these matrices, so that only small matrices are stored and involved in the computation. Numerical experiments support the above model considerations and demonstrate high accuracy in approximating the Classical Scaling embedding results. Using only samples, the approximated embedding results are similar to the ones obtained without any approximation, with an average error of up to .
The paper is organized as follows. In Section 2 we briefly review the Classical Scaling method. In Section 3 we develop two methods for geodesic distances interpolation from a few samples, formulated as a multiplication of smallsized matrices. Section 4 discusses how to select the samples. Next, in Section 5 we reformulate Classical Scaling through the smallsized matrices obtained by the interpolation. Section 6 provides support for the proposed methods with experimental results, by comparing to other methods. Finally, before concluding, Section 7 demonstrates how the proposed methods are applied to two modern applications, namely, texture mapping and weighted averages on surfaces, and discusses an extension of embedding onto a sphere.
2 Classical Scaling
Given points equipped with some similarity measures between them, Multidimensional Scaling (MDS) methods aim at finding an embedding in a low dimensional Euclidean space such that the Euclidean distances are as close as possible to . When the points lie on a manifold, the affinities can be defined as the geodesic distances between them. In this case, which we focus on in this paper, are invariant to isometric deformations of the surface, thus revealing its intrinsic geometry. One way to define this problem, termed Classical Scaling, is through the minimization
(1) 
where , and , and as usual, for and for all . Denoting by the thin eigenvalue decomposition of the symmetric matrix , with only the largest eigenvalues and corresponding eigenvectors, the solution for this problem is given by . This requires the computation of all pairwise distances in the matrix , which is not practical to obtain when dealing with more than a few thousand points. In this paper, we reconstruct from a small set of its columns. We formulate the reconstruction as a multiplication of smallsized matrices, and show how to solve the above minimization problem without explicitly computing the reconstruction. Thus, we overcome the bottleneck of Classical Scaling and significantly reduce both space and time complexities from quadratic to quasilinear.
3 Distances Interpolation
Let be a manifold embedded in some space. Denote by the squared geodesic distance between . In the discrete domain is represented by points, and is represented by a matrix. In this section we develop two methods for reconstructing by decomposing it into smallersized matrices .
3.1 A smoothness based reconstruction
The task of finding smallsized matrices that, when multiplied together, approximate a matrix has been addressed before, see for example [10, 20]. Compared to such classical decomposition, better approximations can be obtained when considering a specific data model. Here, we exploit the fact that the elements of the matrix to be reconstructed are squared geodesic distances, derived from some smooth manifold.
Starting from the continuous case, let be an arbitrary point on the manifold, and denote by the squared geodesic distance from to all other points . Assume the values of are known in a set of samples , and denote these values by . We expect the distance function to be smooth on the manifold, as nearby points should have similar values. Therefore, we aim to find a function that is both smooth and satisfies the constraints . This can be formulated through the minimization
(2) 
where the energy is some smoothness measure of on . One possible, yet somewhat problematic, choice of is the Dirichlet’s energy
(3) 
used in [1], where is the infinitesimal volume element, and the gradient is computed with respect to the manifold. Here, instead, we use a different measure, termed Laplacian energy, that yields a decomposition which is both simple and more accurate, as we will see next. Denoting by the result of the LaplaceBeltrami operator on applied to , this energy is defined as
(4) 
When dealing with two dimentional triangulated surfaces, denote by the vectors and the discretization of and , and define the diagonal matrix such that its diagonal is a discretization of about each correspoding vertex of the triangulations. Denote by the discretization of the LaplaceBeltrami operator, such that . We use the cotangent weights Laplacian in our experiments [22], [3]. However, any discretization matrix of the LaplaceBeltrami can be used. Following these notations, the energy now reads
(5) 
and the interpolation becomes
(6) 
where is a vector holding the values . Recall, is a subset of . Hence, the matrix can be easily defined such that . Since the constraints may be contaminated with noise, a relaxed form of the problem using a penalty function is used instead. This can be formulated as
(7) 
where is a sufficiently large scalar. This is a quadratic equation and its solution is given by
(8)  
(9) 
Recapping, given a set of samples of the distance function , is a reconstruction of in a sense of being a smooth interpolation from its samples. Figures 1 and 2 demonstrate the reconstruction of from its samples for flat and curved manifolds. Recall, is the squared geodesic distance from a point to the rest of the points on the manifold. In this example we chose at the center of the manifold . In Figure 1 the surface is flat, and hence the function is simply the squared Euclidean distance from a point in , , where and are the Euclidean coordinates. We compare the suggested Laplacian energy to the Dirichlet’s energy. As can be seen, when using the Dirichlet’s energy, the function includes sharp discontinuities at the constraint points.
Next, since is actually a column of the matrix , we reconstruct all the other columns in the same manner. Notice that the matrix does not depend on the values of , but only on the samples locations. Hence, it can be computed once and used to reconstruct all columns. Moreover, it is possible to formulate the reconstruction of all the columns of simultaneously by a simple matrix multiplication. Let us choose columns of and let the matrix hold these columns. The locations of the columns in that correspond to the locations of the samples in the manifold. Each column of corresponds to a vector , which contains distances from the samples to a point which corresponds to the column. The columns selection and the computation of are discussed in Section 4. Now, all columns of can be interpolated simultaniously by the multiplication
(10) 
Since is symmetric by definition, we symmetrize its reconstruction by
(11) 
Notice that is a lowrank reconstruction of with a maximal rank of . This is because is of size which means the rank of is at most . Therefore, the rank of is at most . Consequently, can be written as a multiplication of matrices that are not larger than . These matrices can be obtained as follows. Denote by the horizontal concatenation of and , and define the block matrix , where is the identity matrix. It can be verified that
(12) 
Note that we actually only need to keep the matrices and instead of the whole matrix , which reduces the space complexity from to . In our experiments, samples were enough for accurate reconstruction of , where the number of vertices was in the range of to .
In this subsection we derived a lowrank approximation to the matrix from a set of known columns. However, the same interpolation ideas can be used for other tasks, such as matrix completion, in which a matrix is reconstructed from a set of known entries not necessarily arranged in a few columns.
3.2 A learning based reconstruction
In the previous subsection we derived a reconstruction by exploiting the smoothness assumption of the distances on the manifold. Although the laplacian energy minimization results in a simple lowrank reconstruction of , this might not be the best way to interpolate the distances. In this subsection, we formulate a different decomposition based on a learning approach.
As before, we choose arbitrary columns of and let hold these columns. In the previous section, we constructed an interpolation matrix using the LaplaceBeltrami operator as a smoothness measure, such that , or similarly . Here we try to tackle the following question: Is it possible to construct a matrix that would yield a better reconstruction? The matrix for the best reconstruction of in terms of Frobenius norm is obtained by
(13) 
and the solution is given by . In this case, the reconstruction is nothing but a projection of onto the subspace spanned by the columns of ,
(14) 
However, we cannot compute the best since we do not know the entire matrix . Hence, we suggest to learn the coefficients of from the part of that we do know. This can be formulated using a hadamard product as
(15) 
where is a mask of the location of in (the known rows of ). The known rows of can be thought of as our training data for learning . Let hold the rows of that correspond to the known rows of (see Figure 2(a)). The above equation is equivalent to
(16) 
In this form, the number of coefficients in is large with respect to our training data, and hence this will result in overfitting. For instance, for typical with independent columns, the reconstruction of the training data is perfect. Therefore, we next regularize this equation by reducing the number of learned coefficients, in order to avoid overfitting to the training data. Let hold a subset of columns of , and let hold the rows of that correspond to the known rows of . This is demonstrated in Figure 2(b).
The regularized formulation reads
(17) 
where is now a smaller matrix of size , and the solution is given by
(18) 
yielding what is known as a CUR decomposition of ,
(19) 
CUR decompositions are widely used and studied, see for example [10, 20]. Notice that when we obtain
(20) 
which is known as the Nyström decomposition for positive semidefinite matrices [2]. In Figure 4 (CUR) we demonstrate the reconstruction error of with respect to , where . Notice the significant potential improvement over the previous unregularized case of , which results in the Nyström method.
In order to reduce overfitting, we defined the matrix as a subset of the matrix . However, a different choice of can lead to a better reconstruction of . Denote by the eigenvalue decomposition of and denote by the thin eigenvalue decomposition of with only the eigenvalues with largest magnitude and corresponding eigenvectors. The choice of
(21) 
formulates as a linear combination of the columns of instead of a subset of , thus exploiting better learning of the data. Moreover, this will result in a symmetric decomposition of as we will see next. Plugging and in Equation (19) we obtain
(22) 
which can be reduced to
(23) 
Finally, by defining and we obtain the desired decomposition
(24) 
This results in a symmetric lowrank decomposition of with rank (compared to rank in the previous section). It can be thought of as a variant of the Nyström method where is a regularized version of . In Figure 4 (modified) we demonstrate the improvement of this modification. In all our experiments was a good choice for .
4 Choosing the set of columns
The lowrank approximations of developed in Section 3 are constrained to the subspace spanned by chosen columns. Thus, a good choice of columns would be one that captures most of the information of . The farthest point sampling strategy is a 2optimal method for selecting points from a manifold that are far away from each other [14]. The first point is chosen arbitrarily. Then, recursively, the next point is chosen as the farthest (in geodesic sense) from the already chosen ones. After each point is selected, the geodesic distance from it to the rest of the points is computed.
The computation of geodesic distances from a point to the rest can be performed efficiently using the fast marching method [16] for two dimensional triangulated surfaces, and Dijkstra’s shortest path algorithm [8] for higher dimensions. For surfaces with vertices, both methods have complexities of . Therefore, the complexity of choosing points with the farthest point sampling strategy is .
Using the farthest point sampling strategy we obtain samples from the manifold and their distances to the rest of the points, which when squared, correspond to columns of . Since the samples are far from each other, the corresponding columns are expected to be far as well (in sense) and serve as a good basis for the image of . While other columns selection methods need to store the whole matrix or at least scan it several times to decide which columns to choose, here, by exploiting the fact that columns correspond to geodesic distances on a manifold, we do not need to know the entire matrix in advance. The Farthest point sampling method is described in Procedure 1.
5 Accelerating Classical Scaling
In this section, we show how to obtain the solution for Classical Scaling using the decomposition constructed in the previous sections. A straightforward solution would be to plug instead of and repeat the steps: (1) Compute the thin eigenvalue decomposition of with the largest eigenvalues and corresponding eigenvectors. This can be written as
(25) 
(2) Compute the embedding by
(26) 
However, this solution requires storing and computing , resulting in high computational and memory complexities.
To that end, we propose the following alternative. Denote by the thin QR factorization of . Note that the columns of are orthonormal, and is a small upper triangular matrix. We get
(27) 
Denote by the thin eigenvalue decomposition of , keeping only the first eigenvectors and corresponding largest eigenvalues. This can be written as
(28) 
and hence
(29) 
Since is orthonormal as a product of orthonormal matrices and is diagonal, this is actually the thin eigenvalue decomposition of as in Equation (25), obtained without explicitly computing . Finally, the solution is given by
(30) 
We term the acceleration involving the decomposition in Subsection 3.1 as FastMDS (FMDS) and the acceleration involving the decomposition in Subsection 3.2 as NyströmMDS (NMDS), and sum them up in Procedures 2 and 3.
6 Experimental Results
Throughout this section, we refer to the classical scaling procedure (Section 2) as MDS and compare it to its following approximations: The proposed fastMDS (FMDS, Procedure 2), The proposed NyströmMDS (NMDS), (Procedure 3), Spectral MDS (SMDS [1]), LandmarkIsomap [31], Sampled Spectral Distance Embedding (SSDE, [7]), and ISOMAP using Nyström and incremental sampling (ISMDS [33]). All the above methods are mentioned in the introduction. In Figures 8 and 9, we refer by best to the best rank approximation of a matrix with respect to the Frobenious norm. For symmetric matrices this is known to be given by the thin eigenvalue decomposition, keeping only the eigenvalues with largest magnitude and corresponding eigenvectors. Unless specified otherwise, we use the nonrigid shapes from SHREC2015 database [18] for the experiments, and each shape is downsampled to approximately vertices. For SMDS we use eigenvectors. The parameter is set to .
The output of MDS, known as a canonical form [12], is invariant to isometric deformations of its input. We demonstrate this idea in Figure 5, showing nearly isometric deformations of the Giraffe, Paper and Watch shapes and their corresponding canonical forms in obtained using the NMDS method described in Procedure 3. It can be seen that the obtained canonical forms are approximately invariant to the different poses.
Giraffe  

Canonical form  
Paper  
Canonical form  
Watch  
Canonical form 
Unless a surface is an isometric deformation of another flat surface, flattening it into a canonical form would unavoidably involve some deformation of its intrinsic geometry. This deformation is termed as the embedding error, that is defined by the value of the objective function of the MDS,
(31) 
In our next experiment, shown in Figure 6, we compare the canonical forms of the Giraffe and Hand shapes, obtained using the different methods. For presentation purposes we scale each embedding error to and show it below its corresponding canonical form, where is the number of vertices. Qualitative and quantitative comparison show that the proposed methods are similar to MDS in both minimizers and minima of the objective function. Additionally, the embedding results of NMDS are similar to the those of the MDS up to negligible differences.

MDS  NMDS  FMDS  SMDS  Landmark  SSDE 



Stress  5.355  5.355  5.408  5.557  6.288  10.748 


Stress  3.883  3.887  4.047  4.472  9.128  11.191 
In our next experiment, we directly compare the canonical forms of the different methods () to that of MDS (). In order to compare two canonical forms, we align them using ICP [4] and then compute the maximal relative distance among all points, given by , where is the point in the canonical form . Figure 7 shows the results with respect to the number of samples , on the Hand shape. It can be seen that both proposed methods outperform the others in approximating the embedding of MDS.
The proposed methods not only provide fast and accurate approximations to MDS, but can also be used to significantly accelerate the computation of geodesic distances by approximating them. In the next two figures we use the Giraffe shape to show how accurate the approximation is. In Figure 8, we show the average relative reconstruction error of , given by the average of over all , where . Since the compared methods are essentially different lowrank approximations to , we also add the best rank approximation to the comparison.
In Figure 9, we compare a randomly chosen column of to its reconstruction, thus showing the approximation’s distribution over the surface. Figures 8 and 9 show that the proposed methods can be used for accelerating the computation of geodesic distances up to an average error of approximately for FMDS and for NMDS, depending on the number of samples . This can be useful for applications that need to compute and store pairwise geodesic distances.
Finally, in Figure 10, we evaluated the average computation time of MDS, NMDS, FMDS, and SMDS on shapes from the TOSCA database [6], including the computation of the geodesic distances. We change the number of vertices by subsampling and retriangulating the original shapes. The computations were evaluated on a 2.8 GHz i7 Intel computer with 16GB RAM. Due to time and memory limitations, MDS was computed only for small values of .
7 Extensions and applications
In Subsection 7.1, we extend the proposed methods for embedding onto a sphere instead of a Euclidean space. In Subsections 7.2 and 7.3 we demonstrate the proposed methods on two applications. Throughout this section, we use NMDS for computing the embedding.
7.1 Sphere Embedding
So far we discussed MDS as a flattening procedure, meaning that the data is embedded into a flat Euclidean space. Nevertheless, In some cases, it may be beneficial to embed the data onto a sphere, on which distances are measured according to the sphere’s metric.
A dimensional sphere is the set of all points in at an equal distance from a fixed point. Without loss of generality, we assume that this fixed point is the origin, so that the sphere is centered. The geodesic distance between two points on a sphere of radius can be computed by the arc length
(32) 
where is the angle between the vectors in . Given a manifold with pairwise geodesic distances , we aim to find a set of points in , representing the embedding of the manifold on a sphere, such that . Dividing by and applying on both sides leads to
Let be a matrix holding the distances , and define . Since , we can write in matrix formulation
(33) 
Then, we can formulate the problem through the minimization
(34) 
similar to Classical Scaling, where . This minimization can then be solved using the same techniques presented in FMDS or NMDS. Namely, compute the decomposition of with one of the methods in Section 3. Then, follow Section 5 to obtain the truncated eigenvalue decomposition while ignoring the matrix and the scalar . The final solution is then given by .
A step of normalization of the rows of can be added in order to constraint the points to lie on the sphere. Notice that when there exists an exact spherical solution, it will be obtained without this step. This is due to the fact that and therefore . Hence, for an exact solution without embedding errors we get , which is true only when the rows of are normalized.
In the following experiment, a camera was located on a chair at the center of a room, and a video was taken while the chair was spinning. The frames of the video can be thought of as points in a high dimensional space, lying on a manifold with a spherical structured intrinsic geometry. Using the method described in this subsection, we embed this manifold onto a dimensional sphere (without the normalization step), and show the results in Figure 11. A nice visualization of this experiment is added to the supplementary material. As can be seen, both MDS and NMDS result in a similar circular embedding, revealing the intrinsic geometry of the manifold. The pose of the camera for any frame can then be extracted from the embedding, even if the order of the frames in the video is not known. This experiment also demonstrates that the proposed methods are capable of dealing with more complex manifolds in high dimensional spaces and with embedding to metrics other than Euclidean ones.
7.2 Texture Mapping
Texture mapping is a procedure in which flat texture is mapped onto a curved surface. The main problems with most existing methods for texture mapping are that they involve high computational complexities and introduce large distortions to the texture. Zigelman et al. [36] solved this problem using a flat embedding of the surface, obtained with MDS. This embedding defines a mapping of the surface to a plan, on which the texture lies. Unlike other methods, using MDS does not require boundary conditions, and preserves both the local and global structure of the texture. Plugging our acceleration for MDS in this method results in a simple, efficient and accurate texture mapping procedure. In Figure 12 a texture of a chessboard is mapped onto a synthetic surface with positive and negative gaussian curvatures, using both MDS and NMDS. In Figure 13 a texture of stripes and an image of a sheep are projected onto a curved face surface.
7.3 Weighted averages on surfaces
Finding affine combinations in Euclidean space is a fundamental and simple problem. Given anchor points and corresponding weights which sum up to , the affine combination is defined as the weighted average
(35) 
Conversely, given the affine combination, it is required to find the satisfying weights. However, the generalization of weighted averages to surfaces is not trivial, yet if achieved, enables applications such as splines on surfaces, remeshing, shape correspondences, texture transfer and morphing, and more. One way to generalize this problem to surfaces is the Frechet Mean [13], which defines the weighted average as the point that minimizes the sum of weighted squared geodesic distances to the anchors:
(36) 
Panazzo et al. [21] showed that it is possible to get a good estimation for this equation using flat embeddings. According to their method, the surface is first embedded into a Euclidean space using MDS or any other flattening technique. Then, the weighted average is computed as a simple affine combination in the Euclidean space. Since this may result in a point outside of the surface, a third step of projecting it on the surface has to be done. Panazzo et al. combined Metric MultidimensionalScaling with leastsquares meshes [32] for the flattening step. For convenience we term this method MMDS. In MMDS, similar to LandmarkIsomap, a small set of the surface points is chosen and flattened using MDS, and then the rest of the points are interpolated. In the following experiments we analyze the result of plugging the proposed methods in the flattening step.
In figure 14 we compare the canonical forms of the Cat shape from TOSCA database obtained using MDS, NMDS, and MMDS. As for MMDS, it can be seen that problems are expected in extremities such as tail and legs. Using weighted averages on surfaces, geodesic paths between anchors can be computed by correctly adjusting the weights.
In Figure 15 we show the obtained points on the geodesic paths between four selected anchors on the tail of the cat. It can be seen that NMDS significantly improves the results of [21]. In addition, it was enough to use 30 samples for this task in NMDS, where 150 samples used in MMDS were still insufficient, thus reducing complexities by a factor of 5.
Given an initial set of corresponding points on two shapes, weighted averages on surfaces can be used for completing the correspondence for the rest of the points. For each point in the first shape, this is done by computing its weights with respect to the anchors, and then finding the point in the second shape that corresponds to these weights. In figure 16 we compare the correspondence results of Cat0 and Cat2 from TOSCA database, using the protocol from [15]. In both methods the same points were chosen for the initial correspondence, and the same samples were used for the embedding. It can be seen that better results are achieved for NMDS than for MMDS.
8 Conclusions
In this paper we develop two methods for approximating MDS in its Classical Scaling version. In FMDS, following the ideas of spectralMDS, we interpolate the geodesic distances from a few samples using a smoothness assumption. This interpolation results in a multiplication of small matrices that is a lowrank approximation to the distances matrix. Then, time and space complexities of MDS are significantly reduce by reformulating its minimization through the small matrices. The main contribution to the improvement in time complexity and accuracy over SMDS is the fact that we work entirely in the spacial domain, as opposed to SMDS which translates the problem into the spectral domain and truncates the eigenspace. FMDS sets the stage for the second method we term as NMDS. Instead of constructing an interpolation matrix using smoothness assumptions it is now learned from the sampled data, which results in a different lowrank approximation for the distances matrix. A small modification of the method shows its relation to Nyström approximation. Although NMDS performs better than FMDS in all experiments, the interpolation ideas of FMDS can also be used for other tasks, such as Matrix Completion, where instead of being arranged in columns, the known entries might be spread over the matrix. Experimental results show both methods achieve state of the art results in canonical forms computations. In addition, these methods offer an acceleration for computing geodesic distances on manifolds. Two applications show the benefits of being able to accelerate MDS. Finally, An extension to embedding on spheres show these methods can deal with more complex manifold in high dimensional space embedded into nonflat metrics.
9 Acknowledgments
The research leading to these results has received funding from the European Research Council under European Union’s Seventh Framework Programme, ERC Grant agreement no. 267414
References
 [1] Y. Aflalo and R. Kimmel. Spectral multidimensional scaling. Proceedings of the National Academy of Sciences, 110(45):18052–18057, 2013.
 [2] N. Arcolano and P. J. Wolfe. Nyström approximation of wishart matrices. In Acoustics Speech and Signal Processing (ICASSP), 2010 IEEE International Conference on, pages 3606–3609. IEEE, 2010.
 [3] M. Belkin and P. Niyogi. Laplacian eigenmaps and spectral techniques for embedding and clustering. In NIPS, volume 14, pages 585–591, 2001.
 [4] Paul J Besl and Neil D McKay. Method for registration of 3d shapes. In RoboticsDL tentative, pages 586–606. International Society for Optics and Photonics, 1992.
 [5] I. Borg and P. J. Groenen. Modern multidimensional scaling: Theory and applications. Springer, 2005.
 [6] A. M. Bronstein, M. M. Bronstein, and R. Kimmel. Numerical geometry of nonrigid shapes. Springer, 2008.
 [7] A. Civril, M. MagdonIsmail, and E. BocekRivele. SSDE: Fast graph drawing using sampled spectral distance embedding. In Graph Drawing, pages 30–41. Springer, 2007.
 [8] E. W. Dijkstra. A note on two problems in connexion with graphs. Numerische Mathematik l, pages 269–27.
 [9] D. L. Donoho and C. Grimes. Hessian eigenmaps: Locally linear embedding techniques for highdimensional data. Proceedings of the National Academy of Sciences, 100(10):5591–5596, 2003.
 [10] P. Drineas, R. Kannan, and M. W. Mahoney. Fast monte carlo algorithms for matrices iii: Computing a compressed approximate matrix decomposition. SIAM Journal on Computing, 36(1):184–206, 2006.
 [11] H. A. Drury, D. C. Van Essen, C. H. Anderson, C. W. Lee, T. A. Coogan, and J. W. Lewis. Computerized mappings of the cerebral cortex: a multiresolution flattening method and a surfacebased coordinate system. Journal of cognitive neuroscience, 8(1):1–28, 1996.
 [12] A. Elad and R. Kimmel. On bending invariant signatures for surfaces. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 25(10):1285–1295, 2003.
 [13] Maurice Fréchet. Les éléments aléatoires de nature quelconque dans un espace distancié. In Annales de l’institut Henri Poincaré, volume 10, pages 215–310, 1948.
 [14] D. S. Hochbaum and D. B. Shmoys. A best possible heuristic for the kcenter problem. Mathematics of operations research, 10(2):180–184, 1985.
 [15] Vladimir G Kim, Yaron Lipman, and Thomas Funkhouser. Blended intrinsic maps. In ACM Transactions on Graphics (TOG), volume 30, page 79. ACM, 2011.
 [16] R. Kimmel and J. A. Sethian. Computing geodesic paths on manifolds. Proceedings of the National Academy of Sciences, 95(15):8431–8435, 1998.
 [17] T. Kohonen. The selforganizing map. Neurocomputing, 21(1):1–6, 1998.
 [18] Z. Lian, J. Zhang, S. Choi, H. ElNaghy, J. ElSana, T. Furuya, A. Giachetti, R. A. Guler, L. Lai, C. Li, H. Li, F. A. Limberger, R. Martin, R. U. Nakanishi, A. P. Neto, L. G. Nonato, R. Ohbuchi, K. Pevzner, D. Pickup, P. Rosin, A. Sharf, L. Sun, X. Sun, S. Tari, G. Unal, and R. C. Wilson. Nonrigid 3D Shape Retrieval. In I. Pratikakis, M. Spagnuolo, T. Theoharis, L. Van Gool, and R. Veltkamp, editors, Eurographics Workshop on 3D Object Retrieval. The Eurographics Association, 2015.
 [19] R. Liu, V. Jain, and H. Zhang. Subsampling for efficient spectral mesh processing. In Advances in Computer Graphics, pages 172–184. Springer, 2006.
 [20] M. W. Mahoney and P. Drineas. CUR matrix decompositions for improved data analysis. Proceedings of the National Academy of Sciences, 106(3):697–702, 2009.
 [21] Daniele Panozzo, Ilya Baran, Olga Diamanti, and Olga SorkineHornung. Weighted averages on surfaces. ACM Transactions on Graphics (TOG), 32(4):60 2013.
 [22] Ulrich Pinkall and Konrad Polthier. Computing discrete minimal surfaces and their conjugates. Experimental mathematics, 2(1):15–36, 1993.
 [23] R. Pless. Image spaces and video trajectories: Using isomap to explore video sequences. In ICCV, volume 3, pages 1433–1440, 2003.
 [24] S. T. Roweis and L. K. Saul. Nonlinear dimensionality reduction by locally linear embedding. Science, 290(5500):2323–2326, 2000.
 [25] Y. Rubner and C. Tomasi. Perceptual metrics for image database navigation. Springer, 2000.
 [26] B. Schölkopf, A. Smola, and K. Müller. Nonlinear component analysis as a kernel eigenvalue problem. Neural computation, 10(5):1299–1319, 1998.
 [27] E. L. Schwartz, A. Shaw, and E. Wolfson. A numerical solution to the generalized mapmaker’s problem: flattening nonconvex polyhedral surfaces. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 11(9):1005–1008, 1989.
 [28] H. Schweitzer. Template matching approach to content based image indexing by low dimensional euclidean embedding. In Computer Vision, 2001. ICCV 2001. Proceedings. Eighth IEEE International Conference on, volume 2, pages 566–571. IEEE, 2001.
 [29] Gil Shamai, Yonathan Aflalo, Michael Zibulevsky, and Ron Kimmel. Classical scaling revisited. In Proceedings of the IEEE International Conference on Computer Vision, pages 2255–2263, 2015.
 [30] Gil Shamai, Michael Zibulevsky, and Ron Kimmel. Accelerating the computation of canonical forms for 3d nonrigid objects using multidimensional scaling. In Proceedings of the 2015 Eurographics Workshop on 3D Object Retrieval, pages 71–78. Eurographics Association, 2015.
 [31] V. D Silva and J. B. Tenenbaum. Global versus local methods in nonlinear dimensionality reduction. In Advances in neural information processing systems, pages 705–712, 2002.
 [32] Olga Sorkine and Daniel CohenOr. Leastsquares meshes. In Shape Modeling Applications, 2004. Proceedings, pages 191–199. IEEE, 2004.
 [33] H. Yu, X. Zhao, X. Zhang, and Y. Yang. ISOMAP using Nyström method with incremental sampling. Advances in Information Sciences & Service Sciences, 4(12), 2012.
 [34] Kai Yu, Tong Zhang, and Yihong Gong. Nonlinear learning using local coordinate coding. In Advances in neural information processing systems, pages 2223–2231, 2009.
 [35] Yin Zhou and Kenneth E Barner. Locality constrained dictionary learning for nonlinear dimensionality reduction. Signal Processing Letters, IEEE, 20(4):335–338, 2013.
 [36] Gil Zigelman, Ron Kimmel, and Nahum Kiryati. Texture mapping using surface flattening via multidimensional scaling. Visualization and Computer Graphics, IEEE Transactions on, 8(2):198–207, 2002.