EuclideanDistanceMatrices

Utilities for working with matrices of squared Euclidean distances.

• D̃,S = complete_distmat(D, W): Fills in missing entries in an incomplete and noisy squared distance matrix. W is a binary mask indicating available values. (Algorithm 5 from the reference below).
• D̃,E = rankcomplete_distmat(D, W, dim): Same as above, but works on larger matrices and is less accurate. (Algorithm 2 from the reference below).
• P = reconstruct_pointset(D, dim) Takes a squared distance matrix or the SVD of one and reconstructs the set of points embedded in dimension dim that generated D; up to a translation and rotation/reflection. See procrustes for help with aligning the result to a collection of anchors.
• R,t = procrustes(X, Y) Find rotation matrix R and translation vector t such that R*X .+ t ≈ Y
• denoise_distmat(D, dim, p=2) Takes a noisy squared distance matrix and returns a denoised version. p denotes the "norm" used in measuring the error. p=2 assumes that the error is Gaussian, whereas p=1 assumes that the error is large but sparse. The robust factorization comes from TotalLeastSquares.jl.
• posterior Estimate the posterior distribution of locations given both noisy location measurements and distance measurements (not squared), see more details below.

Bayesian estimation of locations

With distance measurements

If both noisy position estimates and noisy distance measurements are available, we can estimate the full Bayesian posterior over positions. To this end, the function psoterior is avialable. We demonstrate how it's used with an example, and start by generating some sythetic data:

using EuclideanDistanceMatrices, Turing
N = 10    # Number of points
σL = 0.1  # Location noise std
σD = 0.01 # Distance noise std (measured in the same unit as positions)

P  = randn(2,N)                       # These are the true locations
Pn = P + σL*randn(size(P))            # Noisy locations
D  = pairwise(Euclidean(), P, dims=2) # True distance matrix (this function exoects distances, not squared distances).
Dn = D + σD*randn(size(D))            # Noisy distance matrix
Dn[diagind(Dn)] .= 0 # The diagonal is always 0

# We select a small number of distances to feed the algorithm, this corresponds to only some distances between points being measured
distances = []
p = 0.5 # probability of including a distance
for i = 1:N
    for j = i+1:N
        rand() < p || continue
        push!(distances, (i,j,Dn[i,j]))
    end
end
@show length(distances)
@show expected_number_of_entries = p*((N^2-N)÷2)

Given the locations P and distances (vector of tuples with indices and distances), we can now estimate the posterior:

part, chain = posterior(
    Pn,
    distances;
    nsamples = 2000,
    sampler = NUTS(),
    σL = σL, # This can also be a vector of std:s for each location, see ?MvNormal for alternatives
    σD = σD  # This can also be a vector of std:s for each location, see ?MvNormal for alternatives
)

The returned object part is a named tuple containing all the internal variables that were sampled. The fields are of type Particles from MonteCarloMeasurements.jl, representing the full posterior distribution of each quantity. The interesting fields are part.P which contains the posterior positions, and part.d which contains the estimated distances. The object chain contains the same information as part, but in the form of a Turing.Chain object.

Note that the number of samples in the posterior will not be the same as the number requested by nsamples since Turing automatically drops bad samples etc.

We can verify that the estimated locations are closer to the true locations than the ones provided by the measurements alone, and plot the results

norm(mean.(part.P) - P) < norm(Pn - P)

scatter(part.P[1,:], part.P[2,:], markersize=6)
scatter!(P[1,:], P[2,:], lab="True positions")
scatter!(Pn[1,:], Pn[2,:], lab="Measured positions")

Under the hood, Turing.jl is used to sample from the posterior. If you have a lot of points, it will take a while to run this function. If the sampling takes too long time, you may try estimating an MAP estimate instead. To do this, run using Optim and then pass sampler = MAP(). More docs on MAP estimation is found here.

With TDOA measurements (time difference of arrival, i.e., differences of distances)

In this setting, we add one location in the matrix of locations, corresponding to the location of the source that generated the ping. We then set the keyword tdoa=true when calling posterior, and let the vector of (i, j, dist) instead be (i,j,tdoa). Below is a similar example to the one above, but adapted to this setting.

N      = 10                      # Number of points
# The standard deviations below can also be supplied as vectors with one element per location
σL     = 0.1                     # Location noise std
σD     = 0.01                    # TDOA noise std (measured in the same unit as positions)
P      = 3randn(2, N)            # These are the true locations
source = randn(2)                # The true source location
Pn     = P + σL * randn(size(P)) # Noisy locations
tdoas  = []
noisy_tdoas = []
p = 0.5 # probability of including a TDOA
for i = 1:N
    for j = i+1:N
        if rand() < p
            di = norm(P[:, i] - source) # Distance from source to i
            dj = norm(P[:, j] - source) # Distance from source to j
            tdoa = di - dj              # This is the predicted TDOA given the posterior locations
            push!(tdoas, (i, j, tdoa))
            push!(noisy_tdoas, (i, j, tdoa + σD * randn()))
        end
    end
end
@show length(tdoas)
@show expected = p * ((N^2 - N) ÷ 2)


part, chain = posterior(
    [Pn source], # We add the source location to the end of this matrix
    noisy_tdoas;
    nsamples = 2000,
    sampler  = NUTS(),
    σL       = σL,
    σD       = σD, # This can also be a vector of std:s for each location, see ?MvNormal for alternatives
    tdoa     = true, # Indicating that we are providing TDOA measurements
)

norm(mean.(part.P[:, 1:end-1]) - P) < norm(Pn - P)

Once again, we visualize the resulting estimate

scatter(part.P[1, 1:end-1], part.P[2, 1:end-1], markersize = 6)
scatter!(P[1, :],  P[2, :],  lab = "True positions")
scatter!(Pn[1, :], Pn[2, :], lab = "Measured positions")
scatter!(
    [part.P[1, end]],
    [part.P[2, end]],
    m = (:x, 8),
    lab = "Est. Source",
)
scatter!([source[1]], [source[2]], m = (:x, 8), lab = "True Source") |> display

Relative vs. Absolute estimates

The function posterior estimates the absolute positions of the sensors in the coordinate system used to provide the location measurements. Oftentimes, the relative positions between the sensors are sufficient, and are also easier to estimate. Estimates of the relative positions are available in the resulting samples from the posterior distribution, but hidden within the samples. If we draw 2000 samples from the posterior, the absolute coordinates of each sample can be aligned to the mean of all samples (using procrustes), after which 2000 samples of the relative positions are available. This relative estimate will have lower variance than the absolute estimate. To facilitate this alignment, we have the function

P_relative = align_to_mean(part.P)
tr(cov(vec(part.P))) > tr(cov(vec(P_relative))) # Test that the covariance matrix is "smaller"

Installation

using Pkg
pkg"add https://github.com/baggepinnen/EuclideanDistanceMatrices.jl"

References

Most of the algorithms implemented in this package are described in the excellent paper "Euclidean Distance Matrices: Essential Theory, Algorithms and Applications" Ivan Dokmanic, Reza Parhizkar, Juri Ranieri and Martin Vetterli
https://arxiv.org/pdf/1502.07541.pdf