
\section{Proofs}

Throughout this section we fix the following notation.
Let $D$ be a distance matrix with squared entries.
Let $r$ be the dimension into which we are embedding. Let $\lambda_1 \leq \hdots \leq \lambda_{n-1}$ be the eigenvalues of $\hat{D}$ and $U$ be the eigenvectors.
Let $\boldsymbol{\lambda}$ be an $n$-dimensional vector where $\boldsymbol{\lambda}_i = \lambda_i\mathbbm{1}_{\lambda_i > 0 \text{ or } i > r}$ for $i=1,\dots,n-1$ and $\boldsymbol{\lambda}_n = 0$.
Let $S = Q*\begin{bmatrix} U & 0 \\ 0 & 1\end{bmatrix}$. 
Let $D_t$ be the solution to Problem \ref{eq:mds} and $D_{\rm{cmds}}$ the resulting EDM from the solution to Problem \ref{eq:cmds}.
Let
\begin{align*}
    C_1 = \sum_{i=1}^{n-1} \boldsymbol{\lambda}^2_i,\ \ \ 
    C_2 = -\sum_{i=1}^{n-1} \boldsymbol{\lambda}_i, \ \ \ 
    C_3 = \frac{n\|(S \circ S)\boldsymbol{\lambda}\|_F^2 - C_2^2}{2}.
\end{align*}

Then the main result of the paper is the following spectral decomposition of the SSTRESS error. 

{\reftheorem{thm:result} If $D$ is a symmetric, hollow matrix then, $\|D_{\rm{cmds}} - D\|_F^2 = C_1 + C_2^2 + C_3$. }

\begin{proof}
First, write $\|D-D_{\rm{cmds}}\|_F^2$ as $\|QDQ-QD_{\rm{cmds}}Q\|_F^2$. Then, the difference between these two matrices can be broken down into three pieces: $\|\hat{D} - \hat{D}_{\rm{cmds}}\|_F^2$ which is given by Lemma \ref{lem:c1} and \ref{lem:conj}, $(\xi(D)-\xi(D_{\rm{cmds}}))^2$ which is given by Lemma \ref{lem:trace}, and 2$\|f(D)-f(D_{\rm{cmds}})\|_F^2$ which is given by Lemma \ref{lem:hollow}. 
\end{proof}



{\reflemma{lem:gram} If $G$ is a positive semi-definite Gram matrix, then $-\frac{1}{2}V \, \text{EDM}(G) \, V = Q\begin{bmatrix} \hat{G} & 0 \\ 0 & 0 \end{bmatrix}Q$}
\begin{proof}
First, note that
\[
    \text{EDM}(G) = \text{diag}(G){\bm{1}}^T -2G + {\bm{1}} \text{diag}(G)^T. 
\]

Then, because $\text{diag}(G){\bm{1}}^T$ is a rank 1 matrix in which every column is the same, when we center it using $V$, we see that $V \text{diag}(G){\bm{1}}^T V = 0$. Similarly, we have that $V{\bm{1}} \text{diag}(G)^T V = 0$. Therefore, 
\[
    V\, \text{EDM}(G) \, V = -2VGV.
\]
Finally, dividing by $-2$ and using the relation between $V$ and $Q$ from Section \ref{sec:matrices}, we obtain the desired result.
\end{proof}

In the following discussion, let $Y_r := {X_{\rm{cmds}}}^TX_{\rm{cmds}}$ and recall that $X_{\rm{cmds}}$ is the solution to the classical MDS problem given in Equation \ref{eq:cmds}.

{\reflemma{lem:c1} 
The value of the objective function obtained by $X_{\rm{cmds}}$ in Equation \ref{eq:cmds} is $\frac{C_1}{4}$. Specifically, we have that
\[  4\|Y_r - (-VDV)/2 \|_F^2 = \sum_{i=1}^{n-1} \boldsymbol{\lambda}^2_i =: C_1.
\]
}
\begin{proof}
Let $B$ be any matrix and note
\begin{align*}
  \left\|B - (-VDV)/2\right\|_F^2 &= \left\|B - Q\begin{bmatrix} -\hat{D}/2 & 0 \\ 0 & 0 \end{bmatrix}Q\right\|_F^2 = \left\|QBQ - \begin{bmatrix} -\hat{D}/2 & 0 \\ 0 & 0 \end{bmatrix}\right\|_F^2.
\end{align*}
The first equality holds because of the relationship between $Q$ and $V$ and the second because of the unitary invariance of the Frobenius norm. In the classical MDS algorithm, we minimize the term on the left hand side with the constraint that $B$ is positive semi-definite with rank at most $r$. Because conjugating a positive semi-definite matrix by unitary matrix results in a positive semi-definite matrix and does not change the rank,  we can solve the MDS problem instead by minimizing the last term with the restriction that $QBQ$ is positive semi-definite with rank at most $r$. We know the solution to this is to keep the $r$ biggest positive eigenvalues of $-\hat{D}/2$. Then, since $\lambda_1 \leq \hdots \leq \lambda_{n-1}$, we have that 
\[
     \|Y_r - (-VDV)/2 \|_F^2 = \frac{1}{4}\sum_{i=1}^{n-1} \boldsymbol{\lambda_i}^2 = \frac{C_1}{4}.
\] 
\end{proof}

{\reflemma{lem:conj}
Suppose $D_{\rm{cmds}} = EDM(Y_r)$, then 
$-\frac{1}{2} \hat{D}_{\rm{cmds}} = \hat{Y_r}$.}
\begin{proof} We see this via the following calculation 
\begin{align*}
    Q\begin{bmatrix} -\hat{D}_{\rm{cmds}}/2 & 0 \\ 0 & 0 \end{bmatrix} Q &= -VD_{\rm{cmds}}V/2 = -V \, \text{EDM}(Y_r) \, V/2 = Q\begin{bmatrix} \hat{Y_r} & 0 \\ 0 & 0 \end{bmatrix}Q.
\end{align*}

The first equality is due to the relationship between $Q$ and $V$ from Section \ref{sec:matrices}. The second is because $D_{\rm{cmds}}$, by definition, is given by $\text{EDM}(Y_r)$. Finally, the last equality is due to Lemma \ref{lem:gram}. 
\end{proof}

{\reflemma{lem:trace}
Suppose $D_{\rm{cmds}} = EDM(X_{\rm{cmds}})$, where $X_{\rm{cmds}}$ is the solution to the cMDS problem given $D$ as input.
If $\text{Tr}(D) = 0$, then 
\[
    (\xi(D) - \xi(D_{\rm{cmds}}))^2 = \left(\sum_{i=1}^{n-1} \boldsymbol{\lambda_i} \right)^2 =: C_2^2.
\]}
\begin{proof}
Since $D_{\rm{cmds}}$ is an EDM, we have that $\text{Tr}(D_{\rm{cmds}}) = 0$ and $\text{Tr}(D-D_{\rm{cmds}}) = 0$. Finally, since $\text{Tr}(QAQ) = \text{Tr}(A)$ for any matrix $A$,
\[
    0 = \text{Tr}(D-D_{\rm{cmds}}) = \text{Tr}(\hat{D} - \hat{D}_{\rm{cmds}}) + \xi(D) - \xi(D_{\rm{cmds}}). 
\]
Then we know from the construction in Lemma 2, that 
\[ 
    \text{Tr}(\hat{D} - \hat{D}_{\rm{cmds}}) =\sum_{i=1}^{n-1} \boldsymbol{\lambda_i}. 
\]
Rearranging and solving gives the desired quantity. 
\end{proof}

{\reflemma{lem:hollow}
\[
    2\|f(D) - f(D_{\rm{cmds}})\|_F^2 =
    \frac{n\|(S\circ S)\boldsymbol{\lambda}\|_F^2 - C_2^2}{2} =: 
    C_3.
\]}
\begin{proof}
First, we note that if $\Lambda$ is the diagonal matrix whose diagonal is given by the first $n-1$ entries of $\boldsymbol{\lambda}$, then 
\[
    U\Lambda U^T = \hat{D} - \hat{D}_{\rm{cmds}}. 
\]
This holds because of Lemmas \ref{lem:gram} and \ref{lem:conj}. Then we can see the following spectral decomposition as well
\[
  Q  \begin{bmatrix}U & 0 \\ 0 & 1 \end{bmatrix}\begin{bmatrix}\Lambda & 0 \\ 0 & 0 \end{bmatrix}\begin{bmatrix}U^T & 0 \\ 0 & 1 \end{bmatrix}Q = Q\begin{bmatrix}\hat{D} - \hat{D}_{\rm{cmds}} & 0 \\ 0 & 0 \end{bmatrix}Q.
\]
\noindent From \citet{HAYDEN1988115}, we have that if $F$ is a symmetric matrix, then there are unique $f$, $\xi$ such that 
\[
    Q\begin{bmatrix} \hat{F} & f \\ f^T & \xi \end{bmatrix}Q
\]
is hollow. The unique $f$, $\xi$ are given by 
\[
    \begin{bmatrix} 2f \\ \xi \end{bmatrix} = \sqrt{n} Q \text{diag}\left(Q \begin{bmatrix} \hat{F} & 0 \\ 0 & 0 \end{bmatrix}  Q \right).
\]

Now, note that $D$ and $D_{\rm{cmds}}$ are already hollow, thus, we have that 
\[
    \begin{bmatrix} 2f(D) \\ \xi(D) \end{bmatrix} = \sqrt{n} Q \text{diag}\left(Q \begin{bmatrix} \hat{D} & 0 \\ 0 & 0 \end{bmatrix} Q \right).
\]
\[
    \begin{bmatrix} 2f(D_{\rm{cmds}}) \\ \xi(D_{\rm{cmds}}) \end{bmatrix} = \sqrt{n} Q \text{diag}\left(Q \begin{bmatrix} \hat{D}_{\rm{cmds}} & 0 \\ 0 & 0 \end{bmatrix} Q \right).
\]
Taking the difference, we get that
\[
    \begin{bmatrix} 2f(D)-2f(D_{\rm{cmds}}) \\ \xi(D) - \xi(D_{\rm{cmds}}) \end{bmatrix} = n Q \text{diag}\!\left(\!Q\! \begin{bmatrix} \hat{D} - \hat{D}_{\rm{cmds}} & 0 \\ 0 & 0 \end{bmatrix} \! Q \! \right).
\]

\noindent Then using the theorem on diagonalization and the Hadamard product from \citet{Million2007TheHP}, we know that 
\[
    \text{diag}\left(Q\begin{bmatrix} \hat{D} - \hat{D}_{\rm{cmds}} & 0 \\ 0 & 0 \end{bmatrix}Q \right) = (S \circ S)\text{diag}(\boldsymbol{\lambda}).
\]
Thus, using Lemma \ref{lem:trace} taking the norm and noting that $Q$ is unitary, we get that
\[
    \left\|\begin{bmatrix} 2f(D)-2f(D_{\rm{cmds}}) \\ \xi(D)-\xi(D_{\rm{cmds}}) \end{bmatrix} \right\|_2^2 = n\| (S \circ S)\text{diag}(\boldsymbol{\lambda}) \|_F^2.
\]
Finally, we have that 
\[
  2\|f(D) - f(D_{\rm{cmds}})\|^2 = \frac{ n\| (S \circ S)\boldsymbol(\boldsymbol{\lambda}) \|_F^2 - C_2^2}{2}.
\]
\end{proof}





Recall $S = Q*\begin{bmatrix} U & 0 \\ 0 & 1\end{bmatrix}$. 

{\reftheorem{thm:reformulate} If $\boldsymbol{c}$ is the solution to Problem \ref{problem:diagonal} and $D_l$ is the solution to problem \ref{problem:relaxed} and if we let $M$ be the $n-1$ by $n-1$ diagonal matrix with the first $n-1$ terms of $\boldsymbol{c}$ on the diagonal, then 
\[ S^TD_lS = \begin{bmatrix} M & f(D) \\ f(D)^T & c_n \end{bmatrix}.
\]}

\begin{proof}

Let $M$ be any matrix and let $Q$ be the unitary matrix in Equation~\ref{eq:Qdef}.  Because the Forbenius norm is invariant under unitary transformations, conjugating both $M$ and $D$ by $Q$ gives us two equivalent measures of the Frobenius distance between two matrices
\[
    \|M - D\|_F = \|QMQ - QDQ\|_F.
\]

\noindent Now we know that $\begin{bmatrix} \hat{D} & f(D) \\ f^T(D) & \xi(D) \end{bmatrix} = QDQ$ and $\hat{D}=U\Lambda U^T$ is the eigen-decomposition of $\hat{D}$. Let us consider the matrix 
\[ 
    R := \begin{bmatrix} U & 0 \\ 0 & 1 \end{bmatrix}.
\]

\noindent Since $R$ is a unitary matrix, we again have that
\[
    \|M - D\|_F = \left\|R^TQMQR - \begin{bmatrix} \Lambda & U^Tf(D) \\ f^T(D)U & \xi(D) \end{bmatrix} \right\|_F.
\]

\noindent Note that $R^TQ = S^T$. Now consider
\[
    QMQ = \begin{bmatrix} \hat{M} & f(M) \\ f^T(M) & \xi(M) \end{bmatrix}.
\]

If we enforce that $M$ is an EDM, then $M$ must be symmetric, hollow, and $\hat{M}$  negative semi-definite. Let us now conjugate by $R$ to obtain 
\[ 
    R^TQMQR = \begin{bmatrix} U^T\hat{M}U & U^Tf(M) \\ f^T(M)U & \xi(M) \end{bmatrix}.
\]
Let $C := U^T\hat{M}U$. Then we have that 
\begin{align*}
     \|M-D\|_F^2 &= \sum_{i\neq j} C_{ij}^2 + \sum_{i=1}^n(C_{ii}-\lambda_i)^2 + 2\|f(M) - f(D)\|_2^2 + (\xi(M)-\xi(D))^2
\end{align*}

\noindent In this case, since $\hat{M}$ is negative semi-definite and $U$ is unitary, we have that $U^T\hat{M}U$ remains negative semi-definite. If we relax the hollow condition on $M$ to $\text{Tr}(M) = 0$, then both $f(M)$ and $\xi(M)$ are unconstrained and we can set $f(M) = f(D)$ and $C_{ij} = 0$ for $i \neq j$. Similarly, we know that we must have $C_{ii} \le 0$ (as $C$ is negative semi-definite).  Therefore, if we relax the condition that $M$ must be hollow, then 

\begin{align*}
    \argmin_{M \in \kappa{r}} \|M-D\|_F^2 &=  \argmin_{\boldsymbol{c} = [c_1, \hdots, c_n]^T, \forall n > i, c_i \le 0, \boldsymbol{c}^T1 = 0, \forall n > j > r c_j = 0} \sum_{i=1}^{n-1}(c_i - \lambda_i)^2 + (c_n-\xi(D))^2.
\end{align*}
\end{proof}

{\reftheorem{thm:lower} If the input matrix $D$ is a symmetric, hollow matrix, then Algorithm \ref{alg:lower} computes $D_l$ the solution to \ref{problem:relaxed}.}

\begin{proof}

To solve the optimization problem given by Problem \ref{problem:relaxed}, it is sufficient to satisfy the KKT conditions. 
To do so, we need to maintain feasibility and find constants $\mu_i > 0$ and $C$ such that $c_i - \lambda_i + \mu_i - C = 0$ for $i = 1, \hdots, r$, for all $n > i > r$ we need $c_i = 0$ and $\xi(M) - \xi(D) - C = 0$. 

First, we note that for $n > i > r$, we have to set $c_i = 0$ which we do in the first for loop. To see that we satisfy the rest, let $C$ be the value of \texttt{sub} on Line 19. Note that, we change \texttt{sub}, when $c_i + \texttt{sub} > 0$. This implies that $c_i > -sub$. Thus, we see that the value of \texttt{sub} only increases. Thus, $C > -C_2/(r+1)$. 

Next we initialize $\mu_i = \lambda_i + C$ if $c_i = 0$ and 0 otherwise. Note if $c_i = 0$ for $i \le r$, then for $j \ge r$, we have that $\lambda_j > 0$. Thus $-C_2$ is positive, thus $C$ is positive. Thus, if $c_i = 0$ for $i \le r$, we have that $\mu_i = \lambda_i + C > 0$. Then Line 13 sets $\xi(M) = \xi(D) + C$. 

If this solution is feasible, we are done. Note now that $\boldsymbol{c}^T1 = -C$. Line 15 checks whether $c_i = \lambda_i + C$ is greater than 0. If it is then this violates a constraint. Hence we set $c_i$ to 0. Doing so may violate the KKT condition $c_i - \lambda_i + \mu_i + C = 0$, so we set $\mu_i$ to be $\lambda_i + C$. This is non negative since $\lambda_i + C \ge \lambda_i + \texttt{sub} \ge 0$. All of these adjustments change the trace so we distribute the excess equally over the rest of the $c_i$. 

Finally, we note that we maintain stationarity by adjusting $\mu_i$ which is initially zero and to which we add only positive amounts. Line 19 maintains trace 0. Thus, we have dual feasibility. Finally, we have complementary slackness as we make $\mu_i \neq 0$ only when we set $c_i = 0$. Thus, $\mu_i c_i = 0$ for all $i$ and we have calculated the optimal solution. 
\end{proof}

{\refprop{prop:bound} If $D$ is a symmetric, hollow matrix, then 
\[
    \|D_l - D\|_F^2 \ge C_1 + \frac{C_2^2}{r+1}
\]}
\begin{proof}

From Theorem \ref{thm:lower}, we have that for if $\lambda_i > 0$ or $i > r$, we set $c_i = 0$.  Recall that 
\[
    C_1 = \sum_{i=1}^{n-1} \lambda_i^2 \mathbbm{1}_{\lambda_i > 0 \text{ or } i > r}
\]

Let $k+1$ be the smallest index we set to 0. Thus, we see that
\[
    \argmin_{\forall n > i, c_i \le 0, \boldsymbol{c}^T1 = 0, \forall n > j > r c_j = 0} \sum_{i=1}^{n-1}(c_i - \lambda_i)^2 + (c_n-\xi(D))^2 = \argmin_{\forall n > i, c_i \le 0, \boldsymbol{c}^T1 = 0} \sum_{i=1}^{k}(c_i - \lambda_i)^2 + (c_n-\xi(D))^2 + C_1
\]

If we now the relax the $c_i \le 0$ constraint for $i = 1, \hdots, n-1$, then the optimal solution is obtained by subtracting $C_2/(r+1)$ from $c_1, \hdots, c_r$, and $c_n$. Thus, we get that
\[
    \|D_l - D\|_F^2 \ge C_1 + \frac{C_2^2}{r+1}.
\]

\end{proof}

{\reflemma{lem:hollow2}
\[
    2\|f(D) - f(D_{\rm{lcmds}})\|_F^2 \le
    \frac{n\|(S\circ S)\boldsymbol{\lambda}\|_F^2 - \frac{C_2^2}{r+1}}{2} 
\]}

\begin{proof}
First, we note that if $\boldsymbol{C}$ is the diagonal matrix whose diagonal is given by the first $n-1$ entries of $\boldsymbol{c}$ then we have that
\[
    U\boldsymbol{C} U^T = \hat{D} - \hat{D}_{\rm{lcmds}}. 
\]
This is holds because of Lemmas \ref{lem:gram} and \ref{lem:conj}. Next, we can see the following spectral decomposition as well
\[
  Q  \begin{bmatrix}U & 0 \\ 0 & 1 \end{bmatrix}\begin{bmatrix}\boldsymbol{C} & 0 \\ 0 & 0 \end{bmatrix}\begin{bmatrix}U^T & 0 \\ 0 & 1 \end{bmatrix}Q = Q\begin{bmatrix}\hat{D} - \hat{D}_{\rm{lcmds}} & 0 \\ 0 & 0 \end{bmatrix}Q.
\]
\noindent From \citet{HAYDEN1988115}, we have that if $F$ is a symmetric matrix, then there are unique $f$, $\xi$ such that 
\[
    Q\begin{bmatrix} \hat{F} & f \\ f^T & \xi \end{bmatrix}Q
\]
is hollow. The unique $f$, $\xi$ are given by 
\[
    \begin{bmatrix} 2f \\ \xi \end{bmatrix} = \sqrt{n} Q \text{diag}\left(Q \begin{bmatrix} \hat{F} & 0 \\ 0 & 0 \end{bmatrix}  Q \right).
\]

Now, note that $D$ and $D_{\rm{lcmds}}$ are already hollow; thus, we have that 
\[
    \begin{bmatrix} 2f(D) \\ \xi(D) \end{bmatrix} = \sqrt{n} Q \text{diag}\left(Q \begin{bmatrix} \hat{D} & 0 \\ 0 & 0 \end{bmatrix} Q \right).
\]
\[
    \begin{bmatrix} 2f(D_{\rm{lcmds}}) \\ \xi(D_{\rm{lcmds}}) \end{bmatrix} = \sqrt{n} Q \text{diag}\left(Q \begin{bmatrix} \hat{D}_{\rm{lcmds}} & 0 \\ 0 & 0 \end{bmatrix} Q \right).
\]
Taking the difference, we get that
\[
    \begin{bmatrix} 2f(D)-2f(D_{\rm{lcmds}}) \\ \xi(D) - \xi(D_{\rm{lcmds}}) \end{bmatrix} = n Q \text{diag}\!\left(\!Q\! \begin{bmatrix} \hat{D} - \hat{D}_{\rm{cmds}} & 0 \\ 0 & 0 \end{bmatrix} \! Q \! \right)
\]

\noindent Next, using the theorem on diagonalization and the Hadamard product from \citet{Million2007TheHP}, we know that 
\[
    \text{diag}\left(Q\begin{bmatrix} \hat{D} - \hat{D}_{\rm{cmds}} & 0 \\ 0 & 0 \end{bmatrix}Q \right) = (S \circ S)\boldsymbol{c}.
\]
Thus, taking the norm and noting that $Q$ is unitary, we get that
\[
    \left\|\begin{bmatrix} 2f(D)-2f(D_{\rm{lcmds}}) \\ \xi(D)-\xi(D_{\rm{lcmds}}) \end{bmatrix} \right\|_2^2 = n\| (S \circ S)\boldsymbol{c} \|_F^2.
\]
Finally, we have that 
\[
  2\|f(D) - f(D_{\rm{lcmds}})\|^2 \le \frac{ n\| (S \circ S)\boldsymbol{c} \|_F^2 - C_2^2/E}{2}.
\]
\end{proof}










\section{Computing the true MDS solution}

The algorithm in \cite{HAYDEN1988115} does not use any dimension constraint and so we use the result from \cite{Qi2014} to adapt it. The algorithm in \cite{HAYDEN1988115} is Bregman's cyclic method or Dykstra's method of alternate projections. Instead of projecting $\hat{D}$ onto the cone of negative semi definite matrices, we project onto $\mathcal{E}(r)$.

We performed each true MDS computation until the change in $\hat{D}$ is at most $0.1$, where the difference is measured as the Frobenius norm. 

\section{Experiments}

First, we note that the experiments do not require extensive computational resources. All experiments were run on a machine with 4 virtual CPUs and 16 GB of RAM. 

\subsection{Perturbed Metrics}

For the perturbed metric, we sampled a Gaussian random matrix $G$, multiplied this by a constant $C$. Then we averaged the noise with its transpose and the main diagonal entries to 0. This noise matrix was then added to the metric $D$ to obtain the input $D_p$. To see the ratio of signal to noise, we looked at $\|D\|_F/\|G\|_F$. For the various datasets, these ratios are as follows. Heart 1.8, MNIST 1.5, FMN!ST 2.1, and CIFAR10 5.8.

Note an alternate way of computing perturbed distances could have been as follows. We add the noise the non-squared distance entries and then square the entries of the new matrix. We ran this experiment as well and we got the following results. These results have the exact same trends as those reported in the main text. 

\begin{figure*}[!h]
    \subfloat[Error decomposition]{\includegraphics[width=0.33\linewidth]{fig-appendix/heart-new-perturbed-err-decomp.png}} 
    \subfloat[Input Error]{\includegraphics[width=0.33\linewidth]{fig-appendix/heart-new-perturbed-input-err-correct.png}}
    \subfloat[Original Error]{\includegraphics[width=0.33\linewidth]{fig-appendix/heart-new-perturbed-orig-err-correct.png}}
    \caption{Plots showing the results on the Heart dataset for the alternate method of perturbation.}
    \label{fig:heart-new}
\end{figure*}

\begin{figure*}[!h]
    \subfloat[Perturbed - MNIST]{\includegraphics[width=0.33\linewidth]{fig-appendix/Perturbed_NN_MNIST.png}} 
    \subfloat[Perturbed - FMNIST]{\includegraphics[width=0.33\linewidth]{fig-appendix/Perturbed_NN_FMNIST.png}}
    \subfloat[Perturbed - CIFAR10]{\includegraphics[width=0.33\linewidth]{fig-appendix/Perturbed_NN_CIFAR.png}}
    \caption{Plots showing the classification accuracy for a 3 layer neural network for the alternate method of perturbation.}
    \label{fig:NN-new}
\end{figure*}



\subsection{Classification}

For classification, we were also interested in the performance of a simpler classifier. So we also tested the 1 nearest neighbor classifier. For this experiment, we were also interested in what true SSTRESS does. Hence, instead of using 2000 data points, we used 500 data points. Then for each of the 500 data points, we classified it using the label of the closest embedded point. The results can be seen in Figure \ref{fig:KNN}

\begin{figure*}[!h]
    \centering
    \subfloat[Missing Data - MNIST]{\includegraphics[width=0.33\linewidth]{fig/Missing_KNN_MNIST.png}} 
    \subfloat[Missing Data - FMNIST]{\includegraphics[width=0.33\linewidth]{fig/Missing_KNN_FMNIST.png}}
    \subfloat[Missing Data - CIFAR10]{\includegraphics[width=0.33\linewidth]{fig/Missing_KNN_CIFAR.png}}\\
    
    \subfloat[Isomap - MNIST]{\includegraphics[width=0.33\linewidth]{fig/Isomap_KNN_MNIST.png}} 
    \subfloat[Isomap - FMNIST]{\includegraphics[width=0.33\linewidth]{fig/Isomap_KNN_FMNIST.png}}
    \subfloat[Isomap - CIFAR10]{\includegraphics[width=0.33\linewidth]{fig/Isomap_KNN_CIFAR.png}}\\
    
    \subfloat[Perturbed - MNIST]{\includegraphics[width=0.33\linewidth]{fig/Perturbed_KNN_MNIST.png}} 
    \subfloat[Perturbed - FMNIST]{\includegraphics[width=0.33\linewidth]{fig/Perturbed_KNN_FMNIST.png}}
    \subfloat[Perturbed - CIFAR10]{\includegraphics[width=0.33\linewidth]{fig/Perturbed_KNN_CIFAR.png}}
    \caption{Plots showing the classification accuracy for the 1 nearest neighbor classifier.}
    \label{fig:KNN}
\end{figure*}



\subsection{Neural Network Details}

We used a 3 layer neural network. The two hidden layers had 100 nodes each. The last layer of the network had a log-softmax layer, and we used the negative log likelihood loss. We trained the networks using ADAM optimizer for 500 epochs. 
\section{Conclusion}
\label{sec:conclusion}
\section{Experiments}
\label{sec:experiments}

In this section, we do two things. First, we empirically verify all of the theoretical claims. Second, we show that on the downstream task of classification, if we use cMDS to embed the data, then as the embedding dimension increases, the classification accuracy gets worse. Thus, suggesting that the embedding quality of cMDs degrades. 

\textbf{Verifying theoretical claims.} To do this, we look at three different types of metrics. First, are metrics that come from graphs and for these we use Celegans \cite{nr} and Portugal \cite{rozemberczki2019multiscale} datasets. Second, are metrics obtained in an intermediate step of Isomap. Third, are perturbations of Euclidean metrics. For both of these metrics we use the heart dataset \cite{Detrano1989InternationalAO}. In all cases, we show that the classical MDS algorithm does not perform as previously expected; instead, it matches our new error analysis. In each case, we demonstrate that our new algorithm Lower + cMDS outperforms cMDS in relation to SSTRESS. 





\textbf{Results.} First, let us see what happens when we perturb an EDM. Here we consider two different measures. Let $D$ be the original input, and  $D_p$ be the perturbed. Then we measure 
\[
   \frac{\|D_p - D_{\rm{cmds}}\|_F^2}{\|D_p\|_F^2}  ~~\text{ and }~~ \frac{\|D-D_{\rm{cmds}}\|_F^2}{\|D\|_F^2}.
\]
As we can see from Figure \ref{fig:input} and \ref{fig:original}, as the embedding dimension increases both quantities eventually increase. The increase in the relative SSTRESS is as we predicted theoretically. The increase in the second quantity suggests that cMDS does not do a good job of denoising either. This suggests that the cMDS objective is not a good one. 

Next, we see how cMDS does on graph datasets and with Isomap on Euclidean data. Here we plot the relative squared error ($\|D-D_{\rm{cmds}}\|_F^2/\|D\|_F^2$).  Figure \ref{fig:final} shows that, as predicted, as the embedding dimension increases, the cMDS error eventually increases as well. For both the Portugal dataset alone and with Isomap on the heart dataset, this error eventually becomes worse than the error when embedding into two dimensions! As we can see from Figure \ref{fig:error}, this increase is \textbf{exactly} due to the $C_2^2$ term. Also, as heuristically predicted, the $C_3$ term is roughly constant.

Let us see how the true SSTRESS solution and new approximation algorithm perform. We look at the relative squared error again. First, Figure \ref{fig:final} shows us that our lower bound tracks the true SSTRESS solution extremely closely. We can also see from Figure \ref{fig:original} that the true SSTRESS solution and the Lower + cMDS solution are closer to the original EDM compared to the perturbed matrix. Thus, they are better at denoising than cMDS. Finally, we see that our new algorithm Lower + cMDS performs better than just cMDS and, in most cases, fixes the issue of the SSTRESS error increasing with dimension. We also see that $\|D_{\rm{lcmds}} - D\|_F^2 - \|D_t - D\|_F^2$ is also roughly constant.










We point out cMDS and Lower+cMDS have comparable running times as they can respectively be computed with one or two spectral decompositions. Computing $D_t$ requires computing a spectral decomposition in every round, with larger datasets requiring hundreds of rounds till they appear to converge. This is a significant advantage of Lower+cMDS over True, and also the reason we were not able to compute True for the Portugal dataset above, and the classifications experiments below.

\textbf{Classification.} For the classification tasks, we switch to more standard benchmark datatsets: MNIST, Fashion MNIST, and CIFAR10. If we treat each pixel as a coordinate then these datasets are Euclidean and we cannot demonstrate the issue with cMDS. We obtain non-Euclidean metrics in three ways. First, we compute the Euclidean metric and then perturb it with a symmetric, hollow matrix whose off diagonal entries are drawn from a Gaussian. Second, we construct a $k$ nearest neighbor graph and then compute the all pair shortest path metric on this graph. This is the metric that Isomap tries to embed using cMDS. Third, we imagine each image has missing pixels and then compute the distance between any pair of images using the pixels that they have in common (as done in \cite{balzano2010column, mrmissing}). Thus, we have 9 different metrics to embed. For each dataset, we constructed all 3 metrics for the first 2000 images. We embedded each of the nine metrics into Euclidean space of dimensions from 1 to 1000 using cMDS and our Lower + cMDS algorithm. We do not embed by solving Equation \ref{eq:mds} as this is computationally expensive. For each embedding, we then took the first 1000 images are training points and trained a feed-forward 3 layer neural network to do classification. Results with a nearest neighbor classifier are in the appendix. We then tested the network on the remaining 1000 points. We can see from Figure \ref{fig:NN}, that as the embedding dimension increases, the classification accuracy drops significantly when trained on the cMDS embeddings. On the other hand, when trained on the Lower + cMDS embeddings, the classification accuracy does not drop, or if it degrades, it degrades significantly less. Thus, showing that the increase in SSTRESS for cMDS on non Euclidean metrics, does result in a degradation of the quality of the embedding and that our new method fixes this issues to a large extent. 

\begin{figure*}[!h]
    \centering
    \subfloat[Missing Data - MNIST]{\includegraphics[width=0.33\linewidth]{fig/Missing_NN_MNIST.png}} 
    \subfloat[Missing Data - FMNIST]{\includegraphics[width=0.33\linewidth]{fig/Missing_NN_FMNIST.png}}
    \subfloat[Missing Data - CIFAR10]{\includegraphics[width=0.33\linewidth]{fig/Missing_NN_CIFAR.png}}\\
    
    \subfloat[Isomap - MNIST]{\includegraphics[width=0.33\linewidth]{fig/Isomap_NN_MNIST.png}} 
    \subfloat[Isomap - FMNIST]{\includegraphics[width=0.33\linewidth]{fig/Isomap_NN_FMNIST.png}}
    \subfloat[Isomap - CIFAR10]{\includegraphics[width=0.33\linewidth]{fig/Isomap_NN_CIFAR.png}}\\
    
    \subfloat[Perturbed - MNIST]{\includegraphics[width=0.33\linewidth]{fig/Perturbed_NN_MNIST.png}} 
    \subfloat[Perturbed - FMNIST]{\includegraphics[width=0.33\linewidth]{fig/Perturbed_NN_FMNIST.png}}
    \subfloat[Perturbed - CIFAR10]{\includegraphics[width=0.33\linewidth]{fig/Perturbed_NN_CIFAR.png}}
    \caption{Plots showing the classification accuracy for a 3 layer neural network.}
    \label{fig:NN}
\end{figure*}



\section{Introduction}
\label{sec:introduction}

\begin{enumerate}
	\item What is MDS?
	\item What's the classical algorithm for solving it? cMDS
	\item What's wrong with MDS as a problem? and cMDS as an algorithm? MDS is difficult to solve, folks have created a zillion variants (by generalizing and changing the objective function). cMDS doesn't actually solve the MDS problem. cMDS as an algorithm is not particularly robust to outliers or perturbations.
	\item EXCEPT no one has analyzed just how non-robust it is! Despite a zillion alternatives, heuristics, and (approximation) algorithms. Cayton and Dasgupta have some heuristic arguments as to how non-robust cMDS and MDS can be in their Robust Euclidean Embeddings paper. They also point out the difference between the optimization problem that cMDS solves and that of MDS.
	\item we seek to quantify this difference -- just how big a problem is this? Are these MDS alternatives and generalizations any better than cMDS? When? What about if we start with data not from a Euclidean space? What if we specify a low-rank or a high-rank?
\end{enumerate}


Learning an embedding into Euclidean space is a standard pre or post processing step in many data pipelines. One of the most common ways this is done, is using the classical Multidimensional Scaling Algorithm (cMDS). For cMDS, given a distance matrix $D_p$, we learn an embedding $X_{cmds}$ whose columns are vectors representing the embedding by solving the following problem. 
\[
    X_{cmds} = \argmin_{X \in \mathbb{R}^{r \times n}} \left\|X^TX - \left(\frac{-V(D\circ D)V}{2}\right)\right\|,
\]
where $V$ is the centering matrix given by $V := I - \frac{1}{n}J$. Here $I$ is the identity matrix and $J$ is the matrix of all ones.  

Now given an embedding $X$, let $EDM(X)$ be the corresponding Euclidean distance matrix. That is, 
\[
    EDM(X)_{ij} = \|X_i - X_j\|^2,
\]
where $X_i,X_j$ are the $i$th and $j$th columns of $X$. Then we see that $D_{cmds} := EDM(X_{cmds})$ is not the solution to the following problem. 
\begin{equation}
\label{eq:mds}    \min_{D = EDM(X), X \in \mathbb{R}^{r \times n}} \|D - (Dp \circ Dp)\|^2.
\end{equation}
The problem in Equation \ref{eq:mds} is the true Multidimensional Scaling problem. Hence we see that the cMDS algorithm is just an approximation algorithm for the MDS problem. While it is possible to solve the true MDS problem exactly (it is a convex optimization problem), it is computational expensive to do so. Thus, the question that we are interested in as follows. Let $D_{t}$ be the solution to the problem in Equation \ref{eq:mds}, then we are interested in understanding the quantity $err := \|D_t - D_{cmds}\|$. In particular, can we predict when $err$ is small without having to compute $D_t$. Doing so provide practitioners with multiple advantages. 

\begin{enumerate}
    \item If $err$ is guaranteed to be small, then we can use the cMDS algorithm without having to worry about loss in quality of the solution. 
    \item If $err$ is big, we can now make an informed decision about whether benefits of the speed of cMDS algorithm versus the quality of the solution. 
    \item Understanding when $err$ is big helps us design algorithms to approximate the solution to MDS problem that perform better when cMDS fails. 
    \item Understanding $err$ is also the first step in rigorously quantifying the robustness of the cMDS algorithm.
\end{enumerate}

\section{Introduction}
\label{sec:introduction}


Multidimensional scaling (MDS) refers to a class of techniques for analyzing pairwise dissimilarities between arbitrary objects. The goal is to find an embedding of points
in a (usually) low dimensional space where the pairwise distances between these points
represent, as best as possible, the given pairwise distances between objects \cite{carroll1998multidimensional, borg2005modern}.
Apart from the
general usefulness of dimensionality reduction, MDS has been used in a wide
variety of applications including data visualization, data preprocessing,
network analysis, bioinformatics,
 and data
exploration.  Due to its long history and being well studied, MDS has many
variations such as nonmetric MDS \cite{shepard1962analysis,
shepard1962analysis2}, multi-way MDS \cite{kroonenberg2008applied}, multi-view
MDS \cite{bai2017multidimensional}, confirmatory or constrained MDS
\cite{heiser1983constrained}, etc. (See \cite{france2010two,
cox2008multidimensional} for surveys).

The basic MDS formulation involves minimizing an objective function over a space
of embeddings. There are two main objective functions associated with MDS: STRESS and STRAIN.
The STRAIN objective (see Equation \ref{eq:cmds} below) was introduced by Torgerson
\cite{torgerson1952}, whose algorithm to solve for this objective is now commonly referred to as the classical MDS algorithm (cMDS).
This algorithm first centers the squares of the given distance matrix and then uses
its spectral decomposition to extract the low dimensional embedding.
cMDS is one of the oldest and most popular methods for MDS, and its popularity is in part due to the fact that this decomposition is very fast and can scale to large 
matrices. The point set produced by cMDS, however, is not necessarily the point set whose Euclidean distance matrix is closest under say Frobenius norm to the input dissimilarity matrix. This type of objective is instead captured by STRESS, which comes in a variety of related forms. In particular, in this paper we consider the SSTRESS objective, shown in Equation \ref{eq:mds} below (and we remark below on its relation to other types of STRESS). There are algorithms to solve the SSTRESS objective, but the main drawback is that they are slow in comparison to cMDS (\cite{takane1977nonmetric,HAYDEN1988115,Qi2014}). Thus, many practitioners default to using cMDS and do not optimize for SSTRESS.


In this paper, we shed light on the theoretical and practical differences
between optimizing for these two objectives. 
Specifically, let $err$ denote the value obtained by applying the SSTRESS objective function to the optimal STRAIN solution. 
Our main theorem, Theorem~\ref{thm:result}, decomposes $err$ into three components.
This decomposition gives insight into when and why cMDS can fail with respect to the SSTRESS objective. In particular, for Euclidean inputs, $err$ naturally decreases as the embedding dimension increases. For non-Euclidean inputs, however, our decomposition shows that after an initial decrease, counterintuitively $err$ can actually increase as the embedding dimension increases. 
This surprising increasing behavior is shown empirically on several real world data sets. 
In practice one may not know a priori what dimension to embed into, though might assume it suffices to embed into some sufficiently large dimension. Importantly, these results demonstrate that when using cMDS to embed into lower dimensions, choosing a dimension too large can actually increase error. Moreover, our decomposition can be used to quickly predict the dimension where this increase might occur. 




This degradation of the cMDS solution is of particular concern in relation to the robustness in the presence of noisy or missing data, as may often be the case for real world data. Several authors \cite{cayton2006robust, mandanas2016robust, forero2012sparsity} have proposed proposed variations to specifically address robustness with cMDS. However, our decomposition of $err$, suggests a novel approach. Specifically, by attempting to directly correct for the problematic term in our decomposition (which resulted in $err$ increasing with dimension) we produce a new lower bound solution. We show empirically that this lower bound corrects for $err$ increasing, both by itself and when used as a seed for cMDS. Crucially the running time of our new approach is comparable to cMDS, rather than the prohibitively expensive optimal SSTRES solution. Finally, and perhaps more importantly, we show that if we add noise or missing entries to real world data sets, then our new solution outperforms cMDS in terms of the downstream task of classification accuracy, under various classifiers.





\textbf{Problem statements and contributions}

Given a dissimilarity matrix \( D \) whose entries are squared, here we consider the embedding producing the closest EDM to $D$ under Frobenius norm. Specifically, given an embedding $X$, $EDM(X)$ is the corresponding Euclidean distance matrix, that is 
\[
    \text{EDM}(X)_{ij} = \|X_i - X_j\|^2_F,
\]
where $X_i,X_j$ are the $i$th and $j$th columns of $X$. Then we are interested in the matrix, 
\begin{equation}
\label{eq:mds}   D_t := \argmin_{D' = EDM(X), X \in \mathbb{R}^{r \times n}} \|D' - D\|^2_F,
\end{equation}
Note that the reason we assume our dissimilarity matrix has squared entries is because the standard EDM definition uses squared entries (see further discussion of EDM's below). 
Equation~\ref{eq:mds} is a well studied objective \citep{takane1977nonmetric,HAYDEN1988115,Qi2014}, and is commonly referred to as SSTRESS. 
There are a number of similarly defined objectives. 
If one considers this objective when the matrix entries are not squared (i.e.\ $\sqrt{D'_{ij}}-\sqrt{D_{ij}}$), then it is referred to as STRESS. 
If one further normalize each entry of the matrix difference by the input distance value (i.e.\ $(\sqrt{D'_{ij}}-\sqrt{D_{ij}})/\sqrt{D_{ij}}$) then it is called Sammon Stress. In this paper, we are less concerned with the differences between different types of Stress, and instead focus on how the cMDS solution behaves generally under a Stress type objective. Thus for simplicity we focus on SSTRESS.

As mentioned above, cMDS does not solve \ref{eq:mds} but instead solves for the following objective
\begin{equation} \label{eq:cmds}
    X_{\rm{cmds}} := \argmin_{X \in \mathbb{R}^{r \times n}} \left\|X^TX - 
            \left(\frac{-V D V}{2}\right)\right\|,
\end{equation}
where $V$ is the centering matrix given by 
\begin{equation}
	V := I - \frac{1}{n}J.
\label{eq:V_def}
\end{equation} 
Here $I$ is the identity matrix and $J$ is the matrix of all ones. This standard MDS objective was introduced in \cite{torgerson1952} and is commonly referred to as STRAIN.



Let $D_{\rm{cmds}} := \text{EDM}(X_{\rm{cmds}})$ and $D_{t}$ be the
solution to the problem in Equation \ref{eq:mds}. We are interested in
understanding the quantity 
\begin{equation}
	\text{err} := \|D - D_{\rm{cmds}}\|.
\label{eq:err_def}
\end{equation}
Doing so will provide practitioners with multiple advantages and will guide the development of better algorithms. In particular, 





\begin{enumerate}[nosep]
	\item Understanding $\text{err}$ is the first step in rigorously quantifying the robustness of the cMDS algorithm.
    \item If $\text{err}$ is guaranteed to be small, then we can use the cMDS algorithm without having to worry about loss in quality of the solution. 
    \item If $\text{err}$ is big, we can make an informed decision about the benefits of the speed of the cMDS algorithm versus the quality of the solution. 
    \item Understanding when $\text{err}$ is big helps us design algorithms to approximate $D_t$ that perform better when cMDS fails. 
    
\end{enumerate}






The main contribution of our paper are as follows. 
\begin{enumerate}[nosep]
	\item We decompose the error in Equation~\ref{eq:err_def} into three terms that depend on the eigenvalues of a matrix obtained from $D$. Using this analysis, we show that there is a term that tells us that as we increase the dimension that we embed into, eventually, the error starts increasing. 
	\item We verify, using classification as a downstream task, that this increase in the error for cMDS results in the degradation of the quality of the embedding, as demonstrated by the classification accuracy decreasing. 
	\item Using this analysis, we provide an efficiently computable algorithm that returns a matrix $D_l$ such that if $D_t$ is the solution to Equation~\ref{eq:mds}, then $\|D_l - D\|_F \le \|D_t-D\|_F$, and empirically we see that $\|D_l-D_t\|_F \le \|D_{\rm{cmds}}-D_t\|_F$. 
	\item While $D_l$ is not metric, when given as input to cMDS instead of $D$, it results in solutions that are empirically better than the cMDS solution. This modified procedure results in a more natural decreasing of the error as the dimension increases and has better classification accuracy. 
\end{enumerate}


\section{Introduction}
\label{sec:introduction}

Multidimensional scaling (MDS) refers to a class of techniques for embedding data into Euclidean space given pairwise dissimilarities \citep{carroll1998multidimensional, borg2005modern}. Apart from the
general usefulness of dimensionality reduction, MDS has been used in a wide variety of applications including data visualization, data preprocessing, network analysis, bioinformatics,
and data
exploration.  Due to its long history and being well studied, MDS has many variations such as non-metric MDS \citep{shepard1962analysis,
shepard1962analysis2}, multi-way MDS \citep{kroonenberg2008applied}, multi-view
MDS \citep{bai2017multidimensional}, confirmatory or constrained MDS
\citep{heiser1983constrained}, etc. (See \cite{france2010two,
cox2008multidimensional} for surveys).

The basic MDS formulation involves minimizing an objective function over a space
of embeddings. There are two main objective functions associated with MDS: STRESS and STRAIN.

The STRAIN objective (Equation \ref{eq:cmds} below) was introduced by \cite{torgerson1952}, whose algorithm to solve for this objective is now commonly referred to as the classical MDS algorithm (cMDS).
\begin{equation} \label{eq:cmds}
    X_{\rm{cmds}} := \argmin_{X \in \mathbb{R}^{r \times n}} \left\|X^TX - 
            \left(\frac{-V D V}{2}\right)\right\|.
\end{equation}
Here $V$ is the centering matrix given by 
$
	V := I - \frac{1}{n}J,
\label{eq:V_def}
$
and $I$ is the identity matrix and $J$ is the matrix of all ones. cMDS first centers the squares of the given distance matrix and then uses
its spectral decomposition to extract the low dimensional embedding.
cMDS is one of the oldest and most popular methods for MDS, and its popularity is in part due to the fact that this decomposition is fast and can scale to large 
matrices. The point set produced by cMDS, however, is not necessarily the point set whose Euclidean distance matrix is closest under say Frobenius norm to the input dissimilarity matrix. This type of objective is instead captured by STRESS, which comes in a variety of related forms. In particular, in this paper we consider the SSTRESS objective (see Equation \ref{eq:mds}). 

Specifically, given an embedding $X$, let $EDM(X)$ be the corresponding Euclidean distance matrix, that is 
$\text{EDM}(X)_{ij} = \|X_i - X_j\|^2_F$,
where $X_i,X_j$ are the $i$th and $j$th columns of $X$. If $D$ is a dissimilarity matrix whose entries are squared, then we are interested in the matrix, 
\begin{equation}
\label{eq:mds}   D_t := \argmin_{D' = EDM(X), X \in \mathbb{R}^{r \times n}} \|D' - D\|^2_F.
\end{equation}
Note that the reason we assume our dissimilarity matrix has squared entries is because the standard EDM characterizations uses squared entries (see further discussion of EDMs below). 
Equation~\ref{eq:mds} is a well studied objective \citep{takane1977nonmetric,HAYDEN1988115,Qi2014}. 

There are a number of similarly defined objectives. 
If one considers this objective when the matrix entries are not squared (i.e.\ $\sqrt{D'_{ij}}-\sqrt{D_{ij}}$), then it is referred to as STRESS. 
If one further normalizes each entry of the matrix difference by the input distance value (i.e.\ $(\sqrt{D'_{ij}}-\sqrt{D_{ij}})/\sqrt{D_{ij}}$) then it is called Sammon Stress. In this paper, we are less concerned with the differences between different types of Stress, and instead focus on how the cMDS solution behaves generally under a Stress type objective. Thus for simplicity we focus on SSTRESS. It is important to note that there are algorithms to solve the SSTRESS objective, but the main drawback is that they are slow in comparison to cMDS \citep{takane1977nonmetric,HAYDEN1988115,Qi2014}. Thus, many practitioners default to using cMDS and do not optimize for SSTRESS.


In this paper, we shed light on the theoretical and practical differences between optimizing for these two objectives. 
Let $D_{\rm{cmds}} := \text{EDM}(X_{\rm{cmds}})$, where $X_{\rm{cmds}}$ is the solution to Equation \ref{eq:cmds}, and let $D_{t}$ be the
solution to Equation \ref{eq:mds}. We are interested in
understanding the quantity 
\begin{equation}
	\text{err} := \|D - D_{\rm{cmds}}\|_F^2.
\label{eq:err_def}
\end{equation}
Doing so will provide practitioners with multiple advantages and will guide the development of better algorithms. In particular, 

\begin{enumerate}[nosep]
	\item Understanding $\text{err}$ is the first step in rigorously quantifying the robustness of the cMDS algorithm.
    \item If $\text{err}$ is guaranteed to be small, then we can use the cMDS algorithm without having to worry about loss in quality of the solution. 
    \item If $\text{err}$ is big, we can make an informed decision about the benefits of the speed of the cMDS algorithm versus the quality of the solution. 
    \item Understanding when $\text{err}$ is big helps us design algorithms to approximate $D_t$ that perform better when cMDS fails. 
    
\end{enumerate}


\textbf{Contributions.} Our main theorem, Theorem~\ref{thm:result}, decomposes $\text{err}$ into three components.
This decomposition gives insight into when and why cMDS can fail with respect to the SSTRESS objective. In particular, for Euclidean inputs, $\text{err}$ naturally decreases as the embedding dimension increases. For non-Euclidean inputs, however, our decomposition shows that after an initial decrease, counterintuitively $\text{err}$ can actually increase as the embedding dimension increases. 
In practice one may not know a priori what dimension to embed into, though one might assume it suffices to embed into some sufficiently large dimension. Importantly, these results demonstrate that when using cMDS to embed, choosing a dimension too large can actually increase error. 

This degradation of the cMDS solution is of particular concern in relation to the robustness in the presence of noisy or missing data, as may often be the case for real world data. Several authors \citep{cayton2006robust, mandanas2016robust, forero2012sparsity} have proposed variations to specifically address robustness with cMDS. However, our decomposition of $\text{err}$, suggests a novel approach. Specifically, by attempting to directly correct for the problematic term in our decomposition (which resulted in $\text{err}$ increasing with dimension) we produce a new lower bound solution. We show empirically that this lower bound corrects for $\text{err}$ increasing, both by itself and when used as a seed for cMDS. Crucially the running time of our new approach is comparable to cMDS, rather than the prohibitively expensive optimal SSTRES solution. Finally, and perhaps more importantly, we show that if we add noise or missing entries to real world data sets, then our new solution outperforms cMDS in terms of the downstream task of classification accuracy, under various classifiers. Moreover, our decomposition can be used to quickly predict the dimension where the increase in $\text{err}$ might occur. 

The main contributions of our paper are as follows. 
\begin{enumerate}[nosep]
	\item We decompose the error in Equation~\ref{eq:err_def} into three terms that depend on the eigenvalues of a matrix obtained from $D$. Using this analysis, we show that there is a term that tells us that as we increase the dimension that we embed into, eventually, the error starts increasing. 
	\item We verify, using classification as a downstream task, that this increase in the error for cMDS results in the degradation of the quality of the embedding, as demonstrated by the classification accuracy decreasing. 
	\item Using this analysis, we provide an efficiently computable algorithm that returns a matrix $D_l$ such that if $D_t$ is the solution to Equation~\ref{eq:mds}, then $\|D_l - D\|_F \le \|D_t-D\|_F$, and empirically we see that $\|D_l-D_t\|_F \le \|D_{\rm{cmds}}-D_t\|_F$. 
	\item While $D_l$ is not metric, when given as input to cMDS instead of $D$, it results in solutions that are empirically better than the cMDS solution. In particular, this modified procedure results in a more natural decreasing of the error as the dimension increases and has better classification accuracy. 
\end{enumerate}
\section*{Acknowledgements}
Work on this paper by Gregory Van Buskirk and Benjamin Raichel was partially supported by NSF CAREER Award 1750780.



\nocite{*}


\section{Preliminaries and Background}
\label{sec:prelim}

In this section, we lay out the preliminary definitions and necessary key structural characterizations of Euclidean distance matrices. 

\textbf{cMDS algorithm.} For completeness, we include the classical multidimensional scaling algorithm in Algorithm~\ref{alg:cmds}.

\begin{algorithm}[!ht]
\caption{Classical Multidimensional Scaling.}
\label{alg:cmds}
	\begin{algorithmic}[1]
	\Function{\textsc{cMDS}}{$D$, $r$}
		\State $X = -V*D*V/2$
		\State Compute $\mu_1 \geq \hdots \geq \mu_{r} > 0$, $U$ as the eigenvalues and eigenvectors of $X$
	   \State return $U * \text{diag}(\sqrt{\mu_1}, \hdots, \sqrt{\mu_r}, 0, \hdots, 0)$
	\EndFunction 
	\end{algorithmic}
\end{algorithm}


\textbf{EDM Matrices}
\label{sec:edm}

\begin{defn}\label{def:edm}
$D \in \mathbb{R}^{n \times n}$ is a Euclidean Distance Matrix (EDM) if and only if there exists a $d \in {\mathbb N}$ such that there are points $x_1, \ldots, x_n \in \mathbb{R}^d$ with
\[
    D_{ij} = \|x_i - x_j\|_2^2. 
\]
Note that unlike other distance matrices, an EDM consists of the squares of the (Euclidean) distances between the points.
\end{defn}

This form permits several important structural characterizations of the cone of EDM matrices.
\begin{enumerate}
	\item \cite{Gower_1985,Schoenberg35remarksto} show that a symmetric matrix $D$ is an EDM if only if 
	\[
    	F := -(I-{\bm{1}} s^T)D(I-s{\bm{1}}^T)
	\]
	is a positive semi-definite matrix for all $s$ such that ${\bm{1}}^Ts = 1$ and $Ds \neq 0$.
	\item \citet{schoenberg1938} showed that $D$ is an EDM if and only if $\rm{exp}(-\lambda D)$ is a PSD matrix with $1$s along the diagonal for all $\lambda > 0$. Note here $\rm{exp}$ is element wise exponentiation of the matrix. 
	\item Another characterization is given by \cite{HAYDEN1988115} in which $D$ is an EDM if and only if $D$ is symmetric, has $0$s on the diagonal (i.e., the matrix is hollow), and $\hat{D}$ is negative semi-definite, where $\hat{D}$ is defined as follows
	\begin{equation}
    \label{eq:edmdef}
    	QDQ = \begin{bmatrix} \hat{D} & f \\ f^T & \xi \end{bmatrix}.
	\end{equation}
	Here $f$ is a vector and 
	\begin{equation}
	\label{eq:Qdef}
	     Q = I - \frac{2}{v^Tv}vv^T \,\,\text{for}\,\, v = [1, \ldots, 1, 1+\sqrt{n}]^T.
	\end{equation} 
	 Note $Q$ is Householder reflector matrix (in particular, it's unitary) and it reflects a vector about the hyperplane span$(v)^\perp$.
\end{enumerate}

The characterization from \cite{HAYDEN1988115} is the main characterization that we shall use. Hence we establish some important notation.
\begin{defn} \label{defn:q-decomp} Given any symmetric matrix $A \in \mathbb{R}^{n \times n}$, let us define $\hat{A} \in R^{n-1 \times n-1}, f(A) \in \mathbb{R}^{n-1}$, and  $\xi(A) \in \mathbb{R}$ as follows
\[
    QAQ = \begin{bmatrix} \hat{A} & f(A) \\ f(A)^T & \xi(A) \end{bmatrix}. 
\]
\end{defn}

In addition to characterizations of the EDM cone, we are also interested in the dimension of the EDM. 
\begin{defn} Given an EDM $D$, the dimensionality of $D$ is the smallest dimension $d$, such that there exist points $x_1, \ldots, x_n \in \mathbb{R}^d$ with $D_{ij}=\|x_i-x_j\|_2^2$. 

Let $\mathcal{E}(r)$ be the set of EDM matrices whose dimensionality is at most $r$.
\end{defn}




Using \cite{HAYDEN1988115}'s characterization of EDMs, \cite{Qi2014} show $D \in \mathcal{E}(r)$ if and only if $D$ is symmetric, hollow (i.e., $0$s along the main diagonal), and $\hat{D}$ in Equation \ref{eq:edmdef} is negative semi-definite with rank at most $r$. 

\textbf{Conjugation matrices: $Q$ and $V$.} 
\label{sec:matrices}
Conjugating distance matrices either by $Q$ (in Equation~\ref{eq:Qdef}) or by the centering matrix $V$ is an important component of understanding both EDMs and the cMDS algorithm. 
We observe that $V$ is also (essentially) a Householder matrix like $Q$. Let $w = [1, \ldots, 1]^T$ and observe that $J = ww^T$ so that $V = I - \frac{1}{w^T w} ww^T$. 



\cite{Qi2014} establishes an important connection between $Q$ and $V$ that we make use of in our analysis in Section~\ref{sec:theory}. Specifically, for any symmetric matrix $A$, we have that 
\begin{equation} \label{eq:connection} 
    VAV = Q \begin{bmatrix} \hat{A} & 0 \\ 0 & 0 \end{bmatrix}Q.
\end{equation}
Here $\hat{A}$ is the matrix given by Definition \ref{defn:q-decomp}. This connection gives a new framework in which we can understand the cMDS algorithm. We know that given $D$, cMDS first computes $-VDV/2$. Using Equation \ref{eq:connection}, this is equal to $-Q\begin{bmatrix} \hat{D} & 0 \\ 0 & 0 \end{bmatrix}Q/2$. Thus, when cMDS computes the spectral decomposition of $-VDV/2$, this is equivalent of computing the spectral decomposition of $\hat{D}$. Then using the characterization of $\mathcal{E}(r)$ from \cite{Qi2014}, setting all but the largest 
$r$ eigenvalues to 0 might seem like the optimal solution. However, this procedure ignores condition that the matrix must be hollow. As we shall show, this results in  sub-optimal solutions. 












\section{Theoretical Results}
\label{sec:theory}

Throughout this section we fix the following notation.
Let $D$ be a distance matrix with squared entries.
Let $r$ be the dimension into which we are embedding. Let $\lambda_1 \leq \hdots \leq \lambda_{n-1}$ be the eigenvalues of $\hat{D}$ and $U$ be the eigenvectors.
Let $\boldsymbol{\lambda}$ be an $n$-dimensional vector where $\boldsymbol{\lambda}_i = \lambda_i\mathbbm{1}_{\lambda_i > 0 \text{ or } i > r}$ for $i=1,\dots,n-1$ and $\boldsymbol{\lambda}_n = 0$.
Let $S = Q*\begin{bmatrix} U & 0 \\ 0 & 1\end{bmatrix}$. 
Let $D_t$ be the solution to Problem \ref{eq:mds} and $D_{\rm{cmds}}$ the resulting EDM from the solution to Problem \ref{eq:cmds}.
Let
\begin{align*}
    C_1 = \sum_{i=1}^{n-1} \boldsymbol{\lambda}^2_i,\ \ \ 
    C_2 = -\sum_{i=1}^{n-1} \boldsymbol{\lambda}_i, \ \ \ 
    C_3 = \frac{n\|(S \circ S)\boldsymbol{\lambda}\|_F^2 - C_2^2}{2}.
\end{align*}

Then the main result of the paper is the following spectral decomposition of the SSTRESS error. 

\begin{thm} \label{thm:result} If $D$ is a symmetric, hollow matrix then, $\|D_{\rm{cmds}} - D\|_F^2 = C_1 + C_2^2 + C_3$. 
\end{thm}

The idea behind the proof of Theorem \ref{thm:result} is to decompose
\begin{align*} \label{eq:decompose}
    \hspace{-1em}\|D_{\rm{cmds}} - D\|_F^2 &= \|QD_{\rm{cmds}}Q-QDQ\|_F^2 \\
    &=4\|\hat{D}/2 - \hat{D}_{\rm{cmds}}/2\|_F^2 + (\xi(D)-\xi(D_{\rm{cmds}}))^2 + 2\|f(D)-f(D_{\rm{cmds}})\|_F^2,
\end{align*}
which follows from Definition \ref{defn:q-decomp}.
We relate each of these three terms to $C_1$, $C_2$, and $C_3$.
The following lemmas will work towards expressing each of these terms separately. In the following discussion, let $Y_r := {X_{\rm{cmds}}}^TX_{\rm{cmds}}$ and recall that $X_{\rm{cmds}}$ is the solution to the classical MDS problem given in Equation \ref{eq:cmds}. The proofs for the following lemmas are in the appendix.

\begin{lemma} \label{lem:gram} If $G$ is a positive semi-definite Gram matrix, then $-\frac{1}{2}V \, \text{EDM}(G) \, V = Q\begin{bmatrix} \hat{G} & 0 \\ 0 & 0 \end{bmatrix}Q$
\end{lemma}

\begin{lemma} \label{lem:c1} 
The value of the objective function obtained by $X_{\rm{cmds}}$ in Equation \ref{eq:cmds} is $\frac{C_1}{4}$. Specifically, we have that
$ \displaystyle 4\|Y_r - (-VDV)/2 \|_F^2 = \sum_{i=1}^{n-1} \boldsymbol{\lambda}^2_i =: C_1.$
\end{lemma}

\begin{lemma} \label{lem:conj}
$-\frac{1}{2} \hat{D}_{\rm{cmds}} = \hat{Y_r}$. \end{lemma}

\begin{lemma} \label{lem:trace}
If $\text{Tr}(D) = 0$, then $\displaystyle (\xi(D) - \xi(D_{\rm{cmds}}))^2 = \left(\sum_{i=1}^{n-1} \boldsymbol{\lambda_i} \right)^2 =: C_2^2.$
\end{lemma}

\begin{lemma} \label{lem:hollow} If $D$ is hollow, then 
$\displaystyle 2\|f(D) - f(D_{\rm{cmds}})\|_F^2 =
    \frac{n\|(S\circ S)\boldsymbol{\lambda}\|_F^2 - C_2^2}{2} =: 
    C_3.$
\end{lemma}

The first term in the error is $C_1$. From the definition, we can see $C_1$ is the sum of the squares of the eigenvalues of $\hat D$ corresponding to eigenvectors that are not used (or are discarded) in the cMDS embedding. As we increase the embedding dimension, we use more eigenvectors. Hence as the embedding dimension increases, we see that $C_1$ monotonically decreases. In the case when $D$ is an EDM, we can use all of the eigenvectors so this term will go to zero. On the other hand, if $D$ is not an EDM and $\hat{D}$ has positive eigenvalues (i.e. negative eigenvalues of $-VDV$), then these are eigenvalues that correspond to eigenvectors that cannot be used. Thus, cMDS will always exhibit this phenomenon for $C_1$ regardless of the input $D$ (for both STRAIN and SSTRESS). 

The second term is $C_2^2$. This term looks similar to $C_1$, but instead of summing the eigenvalues squared, we first sum the eigenvalues and then take the square. This subtle difference has a big impact. First, we note that as $r$ increases, $C_2$ becomes more negative. Suppose that $D$ is not an EDM, (i.e., $\hat{D}$ has positive eigenvalues) and let $K$ be the number of negative eigenvalues of $\hat{D}$. Then, since $D$ has trace 0, when $r=K$, $C_2$ is negative. Hence, $C_2$ decreases and is eventually negative as $r$ increases which implies that as $r$ increases, there is a a certain value of $r$ after which $C_2^2$ \textbf{increases. This results in the quality of the embedding decreasing. As we will see, this term will be the dominant term in the total error.}

While $C_1$ and $C_2$ are simple to understand, $C_3$ is more obtuse. To simplify it, we consider the following. If $\delta$ is an entry of a random $n$ by $n$ unitary matrix, then as $n$ goes to infinity the distribution for $\delta\sqrt{n}$ converges to $\mathcal{N}(0,1)$ and the total variation between the two distributions is bounded above by $8/(n-2)$ \citep{easton,Diaconis1987ADD}. Therefore, we can assume that the variance of an entry of a random $n$ by $n$ orthogonal matrix is about $1/n$. So, heuristically, $S \circ S \approx \frac{1}{n} 11^T$ and
\begin{align*}
    n\|(S\circ S)\vec{\lambda}\|^2_F \approx \frac{n}{n^2}\|11^T\boldsymbol{\lambda}\|^2_F 
    = \frac{n}{n^2}\|1C_2\|^2_F 
    = C_2^2.
\end{align*}
Then since $C_3  = \frac{n\|(S\circ S)\boldsymbol{\lambda}\|_F^2 - C_2^2}{2} \approx \frac{C_2^2-C_2^2}{2}$, we see that, at least heuristically, the overall behavior of $C_3$ is dominated by $C_2$.

\begin{algorithm}[!ht]
\caption{Lower Bound Algorithm.}
\label{alg:lower}
	\begin{algorithmic}[1]
	\Function{\textsc{lower}}{$D$, $r$}
		\State Compute $\hat{D}$, $f(D)$ and $\xi(D)$.
		\State Compute $\lambda_1 \leq \hdots \leq \lambda_{n-1}$, $U$ as the eigenvalues and eigenvectors of $\hat{D}$
		\State Initialize $c_i = \lambda_i$ and $c_n = \xi(D)$
		\State Let $\rm{negative}\_C_2 = 0$.
		\For{$i = 1 \hdots n-1$}
		    \If{$\lambda_i > 0$ or $i > r$}
		        \State $\rm{negative}\_C_2 \mathrel{+}= \lambda_i$
		        \State Set $c_i$ to 0
		    \EndIf
	    \EndFor
	    \State $E := $ number of $c$'s not equal to 0
	    \State $\rm{sub} := negative\_C_2/E$
	    \For{$i = 1 \hdots n-1$}
	        \If{$c_{n-i} != 0$}
	            \If{$c_{n-i} + \rm{sub} \le 0$}
	                \State $c_{n-i} \mathrel{+}= \rm{sub}$, $E \mathrel{-}= 1$, 
                    $\rm{negative}\_C_2 \mathrel{-}= \rm{sub}$
	            \Else
	                \State $E \mathrel{-}= 1$, $\rm{negative}\_C_2 \mathrel{+}= c_{n-1}$, $\rm{sub} = negative\_C_2/E$
	                \State $c_{n-1} = 0$
	            \EndIf
	        \EndIf
	   \EndFor
	   \State $c_n \mathrel{+}= \rm{sub}$, $\rm{negative}\_C_2 \mathrel{-}= \rm{sub}$
	   \State Let $\hat{D} = U * \text{diag}(c_1, \hdots, c_{n-1}) *U^T$
	   \State return $Q * \begin{bmatrix} \hat{D} & f(D) \\ f(D)^T & c_n \end{bmatrix}*Q$
	\EndFunction 
	\end{algorithmic}
\end{algorithm}

From the previous discussion, we see that the $C_2^2$ term in the decomposition is the most vexing due to the term  $(\xi(D)-\xi(D_{\rm{cmds}}))^2$. This term is from the excess in the trace; that is, the result of discarding the eigenvalues changes the value of the trace from 0 to non-zero. From Lemma \ref{lem:conj}, we see that cMDS projects $-\hat{D}/2$ onto the cone of PSD matrices (which do not necessarily have trace 0). To retain the trace 0 condition, we perform the following adaptation. Let $\kappa(r)$ be the space of symmetric, trace 0 matrices $A$, such that $\hat{A}$ is negative semi-definite of rank at most $r$ and project $D$ onto $\kappa(r)$ instead. That is, we seek a solution the following problem. 
\begin{equation}
\label{problem:relaxed}
   D_l := \argmin_{A \in \kappa(r)} \|A - D\|_F^2.
\end{equation}

Clearly this problem is a strict relaxation of the SSTRESS MDS problem and hence we get that 
\[
    \|D-D_t\|_F \ge \|D_l - D\|_F
\]
Because we no longer require the solution to be hollow, conjugating by $S$ permits us to rewrite Problem \ref{problem:relaxed} as: 
\begin{equation}
    \boldsymbol{c} := \argmin_{\substack{\boldsymbol{c} \in \mathbb{R}^n, \boldsymbol{c}^T1 = 0, \\ \forall\, n > i,\ c_i \le 0, \\  \forall\, n > j > r,\ c_j = 0}} \sum_{i=1}^{n-1}(c_i - \lambda_i)^2 + (c_n-\xi(D))^2
    \label{problem:diagonal}
\end{equation}

\begin{thm}
\label{thm:reformulate} If $\boldsymbol{c}$ is the solution to Problem \ref{problem:diagonal} and $D_l$ is the solution to problem \ref{problem:relaxed} and if we let $M$ be the $n-1$ by $n-1$ diagonal matrix with the first $n-1$ terms of $\boldsymbol{c}$ on the diagonal, then $S^TD_lS = \begin{bmatrix} M & f(D) \\ f(D)^T & c_n \end{bmatrix}$. 
\end{thm}

Unlike Problem \ref{eq:mds}, Problem \ref{problem:diagonal} can be solved directly. We do so via Algorithm \ref{alg:lower}. We see that in the first loop (lines 6-9), we set $c_i$ to be 0 if $\lambda_i > 0$ or $i > r$ to ensure the proper constraints on $c_i$s (i.e., the eigenvalues). Doing so, we incur a cost of $C_1$. That is, we sum the squares of those eigenvalues. Now the solution after line 9 does not necessarily satisfy $\boldsymbol{c}^T1 = 0$. In fact, at this stage of the algorithm, we have that $\boldsymbol{c}^T1 = C_2$, and we have $E$ many $c$'s that are non zero that we can adjust. That is, we modify these entries on average by $C_2/E$. Because we use a Frobenius norm to measure error, this incurs a cost of $C_2^2E/E^2 = C_2^2/E$. Finally, we note that the major computational step of the algorithm is the spectral decomposition. Hence it has a comparable run time to cMDS. 

\begin{thm} \label{thm:lower} If $D$ is a symmetric, hollow matrix, then Algorithm \ref{alg:lower} computes $D_l$ the solution to \ref{problem:relaxed} in $O(n^3)$ time.
\end{thm}

\begin{prop} \label{prop:bound} If $D$ is a symmetric, hollow matrix, then 
\[
    \|D_l - D\|_F^2 \ge C_1 + \frac{C_2^2}{r+1}
\]
\end{prop}

It is important to note that $D_l$ is not a metric. However, if we want an EDM and hence an embedding, we can use $D_l$ as the input to cMDS algorithm. Noting that $\hat{D_l}$ is negative semi-definite matrix of rank at most $r$, we see that if $D_{\rm{lcmds}}$ is the metric returned by cMDS on input $D_l$, then we have that $ \|D_l - D_{\rm{lcmds}}\|_F^2 = 2\|f(D_l) - f(D_{\rm{lcmds}})\|_F^2 = 2\|f(D) - f(D_{\rm{lcmds}})\|_F^2$. Similar to Lemma \ref{lem:hollow}, we get the following result. 

\begin{lemma} \label{lem:hollow2}
$\displaystyle  2\|f(D) - f(D_{\rm{lcmds}})\|_F^2 =
    \frac{n\|(S\circ S)\boldsymbol{c}\|_F^2 - \frac{C_2^2}{r+1}}{2} =: C_4.$
\end{lemma}

Using Lemma \ref{lem:hollow2}, we see that
\begin{align*}
    \|D_{\rm{lcmds}} - D\|_F^2 &\le \|D_{\rm{lcmds}} - D_l\|_F^2 + \|D_l - D\|_F^2 \\
                          &\le 2\|f(D_{\rm{lcmds}})-f(D_l)\|_F^2 + \|D_t - D\|_F^2.\\
    \Rightarrow \|D_{\rm{lcmds}} - D\|_F^2 - \|D_t - D\|_F^2 &\le  \frac{n\|(S\circ S)\boldsymbol{c}\|_F^2 - \frac{C_2^2}{r+1}}{2}  = C_4
\end{align*}

On the other hand, we upper bound $ \|D_{\rm{cmds}} - D\|_F^2 - \|D_t - D\|_F^2$ using Proposition \ref{prop:bound}. Taking the difference, we get \[
    \|D_{\rm{cmds}} - D\|_F^2 - \|D_t - D\|_F^2 \le  \frac{r}{r+1}C_2^2 + C_3
\]

Using our heuristics for $C_3$ and $C_4$, we expect $D_{\rm{cmds}}$ to be be much further from $D_t$, than $D_{\rm{lcmds}}$. These claims are experimentally verified in the next section. 

\textbf{Solving Equation \ref{eq:mds}}. Another approach would be to solve Equation \ref{eq:mds} directly. However, solving this problem is much more computationally expensive and works such as \cite{Qi2014} present new algorithms to do so. Algorithm \ref{alg:lower} can be used to compute the solution as well. Here we Dykstra's method of alternating projections as \cite{HAYDEN1988115} do, but instead of only projecting $\hat{D}$ on the cone of negative semi definite matrices, without any rank constraint, we project onto $\kappa(r)$ using Algorithm \ref{alg:lower}. The second, projection is then onto the the space of hollow matrices. Alternating these projections gives us an iterative method to compute $D_t$. 




 