Reconstructing a point set from a Euclidean Distance Matrix

Recovering point locations from Euclidean Distance Matrices has practical applications in ultrasound tomography and other fields that must determine the position of recording devices from acoustic signals. This post explains the solution to a concrete problem, drawing on the explanations in Linear Algebra and Learning from Data and a paper on Euclidean distance geometry.

I want to answer the question: How to reconstruct the locations of the original vectors given the Euclidean distance matrix, ?

We assume the distances in satisfy the triangle inequality. This guarantees a position matrix exists, where has rows for points in -dimensional space. After identifying the point set, you can align it with known locations using Procrustes analysis, also known as the Orthogonal Procrustes problem.

Theory

This procedure, known as classical multidimensional scaling (MDS), comes from this paper. The approach extends to real-world problems like noisy data.

Consider a collection of points in a -dimensional Euclidean space. The squared distance between points and is and is calculated using:

Expanding this norm yields:

The matrix equation for the distance matrix is calculated using:

is the column vector of all ones and is a column vector of the diagonal entries of .

The operator is defined which is equivalent to that operates directly on the Gram matrix . This now gives the equation below:

Let the first point be the origin, then the first column of contains the squared norms of the point vectors,

It is now possible to construct the term and its transpose in the equation to calculate .

The Gram matrix can now be found from like so:

The final stage is to identify the point set using Eigenvalue Decomposition (EVD), for example:

Remember that

Example problem

This problem is Q5 from Problem Set from the Linear Algebra and Learning from Data textbook. Given a Euclidean distance matrix, , find the locations of the points:

Let’s start with some useful imports:

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import distance_matrix
from scipy.spatial.distance import cdist

And now let’s define the matrix :

D = np.array([[0, 9, 25], [9, 0, 16], [25, 16, 0]])

Calculate the Gram matrix. To verify it is positive semidefinite, see J. C. Gower, “Euclidean Distance Geometry,” Math. Sci., vol. 7, pp. 1–14, 1982.

G = -0.5 * (D - np.outer(np.ones(3), D[1, :]) - np.outer(D[:, 1], np.ones(3)))
G

array([[9., -0., -0.],
       [-0., -0., -0.],
       [-0., -0., 16.]])

Now use np.linalg.svd to solve .

Q, Lambda, _ = np.linalg.svd(G)
print(Q)
print(Lambda)

[[0. 1. 0.]
 [0. 0. 1.]
 [1. 0. 0.]]
[16. 9. -0.]

Return the original point set using :

np.power(Lambda, 0.5) * Q.T

array([[0., 0., 0.],
       [4., 0., 0.],
       [0., 3., 0.]])

These coordinates satisfy the original distance constraints. Point 1 is at the origin (0, 0), point 2 at distance 4, and point 3 at distance 3, matching the distances specified in the input matrix .