d3 = require("d3@7")
Plot = require("@observablehq/plot@0.6")
Inputs = require("@observablehq/inputs")
// Original coordinates
X = [
[ 1.14623781, 0.1538715 ],
[ 1.52213386, -0.08244349],
[ 1.14854378, -0.50119446],
[ 0.79205965, 0.49071824],
[-0.94213814, 0.37757528],
[-0.13194707, 1.06805394],
[ 0.03096746, 0.154996 ],
[ 0.50951258, 0.93249084],
[ 1.92567077, 0.15466052],
[ 1.42509235, -0.4438636 ],
[ 0.59488429, -0.20794035],
[ 0.72922735, 0.68034407],
[-1.11931582, -0.04090519],
[ 1.95532585, 0.65242416],
[ 0.20594816, 0.84546606],
[ 0.01612213, -0.24717725],
[-0.74877507, 0.67520749],
[ 0.07458595, 0.31634168],
[ 0.88279825, -0.48225769],
[-0.45268068, 0.93194446]
]
// Sliders — no labels
viewof theta = Inputs.range([0, 2 * Math.PI], { value: 0, step: 0.01 })
viewof mu = Inputs.range([-2, 2], { value: 0, step: 0.1 })
// Transform points
transformed = {
const cosT = Math.cos(theta);
const sinT = Math.sin(theta);
return X.map(([x0, y0]) => {
const xRot = cosT * x0 - sinT * y0;
const yRot = sinT * x0 + cosT * y0;
return { x: xRot + mu, y: yRot + mu };
});
}
// Reference points as "x", transformed as circles
Plot.plot({
width: 690,
height: 450,
x: { domain: [-1.5, 2.5] },
y: { domain: [-1, 1.5] },
marks: [
Plot.frame(),
// Larger red "x"
Plot.text(
X.map(([x, y]) => ({ x, y, label: "x" })),
{
x: "x",
y: "y",
text: "label",
fill: "red",
fontSize: 20, // increase size
fontWeight: "bold"
}
),
// Moving configuration
Plot.dot(transformed, { x: "x", y: "y", r: 4, fill: "blue" })
]
})Multidimensional Scaling From Noisy Data
Statistical Learning & Uncertainty Quantification
Siddharth Vishwanath
Based on joint work with
Ery Arias-Castro
Pairwise Dissimilarity Data
Many problems in psychology, sociology, ecology, wireless communication, neuroscience, bioinformatics, etc.
involve relational data in the form of pairwise dissimilarities between items
Items: \(\{1, 2, \ldots, n\}\)
\[
\Delta= (\delta_{ij})_{1 \leq i,j \leq n} \in \mathbb{R}^{n \times n}
\]
Multidimensional Scaling (MDS)
- Given pairwise dissimilarity data \(\Delta= (\delta_{ij}) \in \mathbb{R}^{n \times n}\) and an embedding dimension* \(p \ll n\)
- Find a configuration of points \(\{ \widehat{x}_1, \dots, \widehat{x}_n\} \in \mathbb{R}^p\) such that:
\(\displaystyle\| \widehat{x}_i - \widehat{x}_j \|^2 \approx \delta_{ij} \quad \text{ for all } 1 \leq i,j \leq n\)
*usually \(p=2\) or \(p=3\)
\(\longrightarrow\)
\(\longrightarrow\)
Multidimensional Scaling (MDS)
Classical MDS (CMDS) Algorithm
Given: dissimilarity matrix \(\Delta= (\delta_{ij}) \in \mathbb{R}^{n \times n}\) and \(p \ll n\)
- Double centering:
\(\Delta_c= -\frac{1}{2}H \Delta H\quad\) where \(\quad H = \Big(I_n - \frac{1}{n}\boldsymbol{1}_n \boldsymbol{1}_n^{\mathsf{T}}\Big)\)
- Rank-p spectral decomposition:
\(\Delta_c= {\widehat{U}}\widehat{\Lambda }{\widehat{U}}^{\mathsf{T}}+ \widehat{V}_{\perp} \widehat{\Lambda}_{\perp} \widehat{V}_{\perp}^{\mathsf{T}}\)
- Spectral embedding:
\(\widehat{X}= {\widehat{U}}\widehat{\Lambda}^{1/2}\)
Output: embedding \(\widehat{X}\in \mathbb{R}^{n \times p}\)
\(\Delta\in \mathbb{R}^{n \times n}\)
\(\Delta_c\in \mathbb{R}^{n \times n}\)
\({\widehat{U}}\in \mathbb{R}^{n \times p}\)
\(\widehat{\Lambda}\in \mathbb{R}^{p \times p}\)
\(\widehat{X}\in \mathbb{R}^{n \times p}\)
Dates back to (Young and Householder 1938). Later rediscovered by (Torgerson 1952) and (Gower 1966).
Identifiability of MDS Embeddings
In the realizable setting, there exists a latent configuration \(\color{green}{x_1, \dots, x_n} \in \mathbb{R}^p\) such that \(\color{green}{\delta_{ij} = \|x_i - x_j\|^2}\)
Let \(\color{green}{X} \in \mathbb{R}^{n \times p}\) be the configuration matrix
For any rigid transformation \(g \in \mathcal{G}\) (i.e., rotation + reflection + translation), \(\color{green}{\Delta\big(\color{green}{X}\big)} \equiv \color{green}{\Delta\big(}\color{magenta}{g}(\color{green}{X)\big)}\)
Theorem (Young and Householder 1938; Schoenberg 1935).
In the realizable setting, the CMDS algorithm exactly recovers the latent configuration up to rigid transformations, i.e.,
\[ \mathsf{CMDS}\big(\; \color{green}{\color{green}{\Delta(X)}}, p\;\big) = \color{magenta}{g}(\color{green}{X}) \quad \text{ for some } \color{magenta}{g} \in \mathcal{G} \]
(a short proof is in the appendix)
MDS under noise
In practice, the dissimilarity data is often subject to noise and/or measurement errors, i.e., we observe: \[ \begin{aligned} d_{ij} &= \color{green}{\delta_{ij}} + \color{red}{\varepsilon_{ij}}\\ D &= \underbrace{\color{green}{\Delta}}_{\color{green}{\text{signal}}} + \underbrace{\color{red}{\mathcal{E}}}_{\color{red} {\text{noise}}} \end{aligned} \]
\(=\)
\(+\)
MDS under noise
In practice, the dissimilarity data is often subject to noise and/or measurement errors, i.e., we observe: \[ \begin{aligned} d_{ij} &= \color{green}{\delta_{ij}} + \color{red}{\varepsilon_{ij}}\\ D &= \underbrace{\color{green}{\Delta}}_{\color{green}{\text{signal}}} + \underbrace{\color{red}{\mathcal{E}}}_{\color{red} {\text{noise}}} \end{aligned} \]
- \(\color{green}{X}\) is the configuration matrix underlying the signal \(\Delta\)
- \(\widehat{X}= \mathsf{CMDS}(D, p)\) is the CMDS embedding of \(D\)
The noisy CMDS embedding \(\widehat{X}\) reflects the variability in \(D\) due to \(\color{red}{\mathcal{E}}\)
Example - 1: Sensor Localization
Suppose \(n\) sensors are placed in various locations across the United States, and \[ d_{ij} = \color{red}{\text{noisy}} \text{ squared Euclidean distance between sensors $i$ and $j$} \]
Are the sensors placed in the following locations?
- San Diego, CA
- Phoenix, AZ
- \(\dots\)
- Washington, DC
\[ \mathbb{P}\Big( \text{locations}(d_{ij}) \simeq \{\text{San Diego}, \dots, \text{Washington DC}\} \Big)? \]
Desiderata - II (Uncertainty Quantification)
- For any \(\color{DarkOrange}{\boldsymbol{\alpha}} \in (0, 1)\)
- Can we construct \(\color{DarkOrange}{\boldsymbol{100 \times (1-\alpha)\%}}\) uniform confidence sets \[ \mathcal{C}_{\alpha} = \prod_{i=1}^n \mathcal{C}_{\alpha, i} \subset \mathbb{R}^p \] which provide simultaneous coverage* of the true configuration: \[ \boxed{\phantom{\Bigg()} \mathbb{P}\Big( \exists g \in \mathcal{G}: \; g(\color{green}{x_i}) \in \mathcal{C}_{\alpha, i} \;\;\; \color{Purple}{\forall i \in [n]} \Big) \approx \color{DarkOrange}{\boldsymbol{1 - \alpha}} \phantom{\Bigg()}} \]
*modulo rigid transformations
Illustration
Desiderata - II (Uncertainty Quantification)
- For any \(\color{DarkOrange}{\boldsymbol{\alpha}} \in (0, 1)\)
- Can we construct \(\color{DarkOrange}{\boldsymbol{100 \times (1-\alpha)\%}}\) uniform confidence sets \[ \mathcal{C}_{\alpha} = \prod_{i=1}^n \mathcal{C}_{\alpha, i} \subset \mathbb{R}^p \] which provide simultaneous coverage* of the true configuration: \[ \boxed{\phantom{\Bigg()} \mathbb{P}\Big( \exists g \in \mathcal{G}: \; g(\color{green}{x_i}) \in \mathcal{C}_{\alpha, i} \;\;\; \color{Purple}{\forall i \in [n]} \Big) \approx \color{DarkOrange}{\boldsymbol{1 - \alpha}} \phantom{\Bigg()}} \]
*modulo rigid transformations
Extreme value distributional approximation
\(\frac{T_n - b_n}{a_n}\quad\) vs. \(\quad G\)
\(O\left(\tfrac{\log\log{n}}{\log{n}}\right)\) convergence rate is extremely slow!
Bootstrap confidence sets
Given: \(D = \color{green}{\Delta}+ \color{red}{\mathcal{E}}\)
1. \(\widehat{X}= \mathsf{CMDS}(D, p)\)
2. Compute the residuals \(\color{purple}{E} = D - \Delta(\widehat{X})\)
3. For \(b = 1, \ldots, B\):
3.1 Externally randomize \(\color{purple}{E}\) to get \(\color{red}{\mathcal{E}}^{\flat}\)
• Empirical bootstrap: \(\color{red}{\varepsilon_{ij}} \sim_{\text{iid}} (e_{ij})\) with replacement
• Multiplier bootstrap: \(\color{red}{\varepsilon_{ij}} = e_{ij} \cdot r_{ij}\) where \(r_{ij} \sim_{\text{iid}} N(0, 1)\)
3.2 Generate noisy dissimilarities \(D^{\flat}= \Delta(\widehat{X}) + \color{red}{\mathcal{E}}^{\flat}\)
3.3 Compute \(\widehat{X}^{\flat}= \mathsf{CMDS}(D^{\flat}, p)\)
3.4 Compute \(\widehat T_n^{\flat}= \max_{i}\|\widehat{\Omega}^{-1/2}(\widehat{x}_i - \color{magenta}{\widehat{g}}^{\flat}(\widehat{x}_i^{\flat}))\|\) via Procrustes
4. Compute the bootstrap quantile \(q_{1-\alpha}^{\flat}\) of \(\{\widehat T_n^{\flat}(1), \dots, \widehat T_n^{\flat}(B)\}\)
Return: \(\displaystyle\mathcal{C}_{\alpha, i} := \mathscr{E}\big( \widehat{x}_i,\; \widehat{\Omega}_i,\;\; q^{\flat}_{1-\alpha} \big)\) for \(i=1, \dots, n\)
\(D \in \mathbb{R}^{n \times n}\)
\(\widehat{X}\in \mathbb{R}^{n \times p}\)
\(E \in \mathbb{R}^{n \times n}\)
\(\mathcal{E}^{\flat}\in \mathbb{R}^{n \times n}\)
\(D^{\flat}\in \mathbb{R}^{n \times n}\)
\(\widehat{X}^{\flat}\in \mathbb{R}^{n \times p}\)
Bootstrap distributional approximation
\({T_n}\quad\) vs. \(\quad T_n^{\flat}\)
The bootstrap distributional approximation is substantially better!
Bootstrap Example
\(X\in \mathbb{X}\) is chosen inside an ellipsoid with \(\kappa=2\) and \(d_{ij} = \color{green}{\delta_{ij}} + \color{red}{\varepsilon_{ij}}\) multiplicative/additive noise
\(\color{red}{\varepsilon_{ij}} \sim N(0, 1)\)
\(\color{red}{\varepsilon_{ij}} \sim N(0, \delta_{ij}^2)\)
Why the bootstrap approximation is better
Population world
We observe \(D = \color{green}{\color{green}{\Delta(X)}} + \color{red}{\mathcal{E}}\phantom{{}^{\flat}}\) and \(\widehat{X}= \mathsf{CMDS}(D, p)\phantom{\widehat T_n^{\flat}}\)
Bootstrap world
We observe \(D^{\flat}= \color{green}{\Delta(\widehat{X})} + \color{red}{\mathcal{E}}^{\flat}\) and \(\widehat{X}^{\flat}= \mathsf{CMDS}(D^{\flat}, p)\)
Adaptivity of Multiplier Bootstrap to Heteroscedasticity
\(X\in \mathbb{X}\) is uniformly spaced on a grid and \(d_{ij} = \color{green}{\delta_{ij}} + \color{red}{\varepsilon_{ij}}\) with radial sum/difference noise:
\(\color{red}{\varepsilon_{ij}} \sim N(0, \|x_i\|^2 + \|x_j\|^2)\)
\(\color{red}{\varepsilon_{ij}} \sim N(0, |\|x_i\|^2 - \|x_j\|^2|)\)
Appendix: Comparison
- \(\color{green}{X}\in \mathbb{R}^{n \times p}\) are chosen from an ellipsoid with condition number \(\kappa^2\)
- \(D = \color{green}{\Delta(X)}+ \color{red}{\mathcal{E}}\) under different noise models with \(\color{red}{\varepsilon_{ij}} \sim t_q\) entries
Appendix: Comparison
