Project new point into MDS space

Question

I am trying to project a new point A(x, y, z) into a predefined MDS space in R. This is what I have so far:

set.seed(1)
x <- matrix(rnorm(3*10), ncol = 3)
DM <- dist(x) 
MDS <- cmdscale(DM)

# New data point to be projected
A <- c(1, 2, 3)

I am not including A directly into x then fitting the MDS because it would affect the space coordinates. Is there a practical solution?

---

EDIT
I believe I found a solution by estimating the betas to predict the MDS axis:

x1 <- cbind(1, x) # add intercept
B <- solve(t(x1) %*% x1) %*% t(x1) %*% MDS # Betas

> MDS
             [,1]        [,2]
 [1,] -1.80789362  0.06801597
 [2,] -0.64418055 -0.21163109
 [3,]  0.04694820 -1.27040928
 [4,]  3.39617277 -0.21657115
 [5,] -0.96981358  0.46269025
 [6,] -0.24716695 -0.79861234
 [7,]  0.33620625  0.02618564
 [8,]  0.62473570  1.35544267
 [9,]  0.01895042  0.80023822
[10,] -0.75395865 -0.21534889

> x1 %*% B # same as MDS
             [,1]        [,2]
 [1,] -1.80789362  0.06801597
 [2,] -0.64418055 -0.21163109
 [3,]  0.04694820 -1.27040928
 [4,]  3.39617277 -0.21657115
 [5,] -0.96981358  0.46269025
 [6,] -0.24716695 -0.79861234
 [7,]  0.33620625  0.02618564
 [8,]  0.62473570  1.35544267
 [9,]  0.01895042  0.80023822
[10,] -0.75395865 -0.21534889

A <- c(1, 2, 3)
A <- c(1, A) # add intercept

> A %*% B # coordinates of A in the MDS plane
          [,1]      [,2]
[1,] -2.759456 0.5927178

Is my procedure correct?

Answer 1

If you're using Euclidean distances, then classical MDS is equivalent to PCA, which readily defines a mapping into the low dimensional space, as amoeba mentioned. There should be various threads on this site describing how to do this. Otherwise, De Silva and Tenenbaum (2004) describe how to perform this mapping for classical MDS with arbitrary distances (note that it won't work for non-classical variants of MDS, e.g. non-metric MDS, variants that minimize the stress criterion, etc.). They call this procedure "distance-based triangulation". Although not originally formulated as such, it turns out to work by using the Nyström approximation, which is a way to approximate the eigenvalues/eigenvectors of a large matrix using a smaller submatrix (see Platt 2005).

Suppose we have $n$ training points. The squared distance between points $i$ and $j$ is stored in the $(i,j)$th entry of matrix $\Delta_n$. We use these distances with classical MDS to compute a $k$-dimensional embedding. Let each column of $k \times n$ matrix $L_k$ contain the low-dimensional embedding coordinates of a training point. Let $L_k^\#$ denote the transpose of the pseudoinverse of $L_k$. Instead of starting from scratch, this can be computed using components that were originally used to compute $L_k$ (see the paper for details). Let $\vec{\delta}_i$ denote the $i$th column of $\Delta_n$ (containing the squared distances from point $i$ to all other points), and let $\vec{\delta}_\mu = \frac{1}{n} \sum_{i=1}^n \vec{\delta}_i$ denote the mean of the columns.

Now, suppose we want to map a new point $a$ into the low dimensional space. Compute vector $\vec{\delta}_a$, containing the squared distance from $a$ to every training point. The low dimensional emedding coordinates of $a$ are then given by:

$$\vec{x}_a = -\frac{1}{2} L_k^\# (\vec{\delta}_a - \vec{\delta}_\mu)$$

References:

De Silva and Tenenbaum (2004). Sparse multidimensional scaling using landmark points.

Platt (2005). FastMap, MetricMap, and Landmark MDS are all Nystrom Algorithms.

Answer 2

As per amoeba's comment, another practical solution in R would be to use the stats::prcomp() and its associated predict() method if using the Euclidean distance for your MDS:

set.seed(1)
x <- matrix(rnorm(3*10), ncol = 3)

# New data point to be projected
A <- c(1, 2, 3)

> predict(prcomp(x, rank. = 2), matrix(A, nrow = 1))
           PC1        PC2
[1,] -2.759456 -0.5927178