Mismatched Estimation
in the Distance Geometry Problem

The distance geometry problem (DGP) [1, 2, 3] aims to determine the locations of points in Euclidean space given measurements of the distances between these points. The DGP arises in many important applications. In wireless networks with moving devices in two-dimensional (2D) planes, the positions of these devices can be estimated by solving the DGP [4]. In computational biology, solving the DGP can allow us to determine the locations of atoms in molecules [5]. In robotics, a simple robotic movement can be modeled as a pair of rigid bars rotating along a joint. Whether these bars can flex along a joint to reach a specific point in 3D space can be determined by solving the DGP [6].

The main challenge for the DGP is that pairwise distance measurements are noisy [1]. Because of the noise, it is common to find a set of locations of points that minimizes the sum of squared errors (SSE) between the pairwise distances of a target structure and noisy measurements of these pairwise distances [7]. When the noisy measurements of pairwise distances are Gaussian distributed, minimizing the SSE function yields a maximum likelihood (ML) estimate of a structure of points. However, in practice, it is unlikely that the measurements are exactly Gaussian distributed, and so the SSE function is unlikely to yield an ML estimate. In this work, we will provide better estimation quality by performing ML estimation using the correct noise distribution. We call our approach matched DGP estimation. In contrast, much of the prior art performed mismatched DGP estimation (details below).

Problem formulation. Let $\mathcal{X}=\{\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{N}\}$ be an unobserved set of locations of $N$ points in $K$-dimensional Euclidean space $\mathbb{R}^{K}$. Additionally, let $\mathcal{E}(\mathcal{X})=\{(i,j)\mid(i,j)\in\{1,\dots,N\}^{2},i<j\}$ be the set of edges connecting every pair of points in $\mathcal{X}$. Note that $|\mathcal{E}(\mathcal{X})|=\frac{1}{2}N(N-1)$, where $|\cdot|$ denotes the cardinality of a set. We define the length $d_{ij}^{\mathcal{X}}$ of an edge $(i,j)\in\mathcal{E}(\mathcal{X})$ as the $\ell_{2}$ norm of $\mathbf{x}_{i}-\mathbf{x}_{j}$, $d_{ij}^{\mathcal{X}}=\left\|\mathbf{x}_{i}-\mathbf{x}_{j}\right\|_{2}$. We define $M_{ij}\in\mathbb{N}$, where $\mathbb{N}$ is the set of nonnegative integers, as the number of noisy measurements of the length of edge $(i,j)\in\mathcal{E}(\mathcal{X})$. An example structure $\mathcal{X}$ is shown in Fig. 1. In the DGP, our objective is to determine $\mathcal{X}$ given noisy measurements of lengths of edges in $\mathcal{E}(\mathcal{X})$.

Estimating the locations of the points. We estimate $\mathcal{X}$ using the method of ML. Let the random variable $Y_{ij}^{[m]}\in\mathbb{R}^{+}$ represent the $m$th noisy measurement of the length of edge $(i,j)$, where a realization of $Y_{ij}^{[m]}$ is denoted by $y_{ij}^{[m]}$. For example, $y_{1,2}^{[8]}$ is the $8$th measurement of the length of edge $(1,2)$. We assume that $M_{ij}=M$ for every $(i,j)\in\mathcal{E}(\mathcal{X})$. Let the set of all $Y_{ij}^{[m]}$ be, $\mathscr{Y}=\bigcup_{(i,j)\in\mathcal{E}(\mathcal{X})}\mathscr{Y}_{ij}$, where $\mathscr{Y}_{ij}=\left\{Y_{ij}^{[m]}\right\}_{m=1}^{M}$. The sets $\mathscr{y}$ and $\mathscr{y}_{ij}$ are defined similarly for $y_{ij}^{[m]}$. Our objective is to derive an expression for the likelihood function $f_{\mathscr{Y}}(\mathscr{y}\mid d_{\mathcal{E}(\mathcal{X})})$, where $d_{\mathcal{E}(\mathcal{X})}=\left\{d_{ij}^{\mathcal{X}}\right\}_{(i,j)\in\mathcal{E}(\mathcal{X})}$.

Suppose that the measurements for edge $(i_{1},j_{1})$ are independent of the measurements for edge $(i_{2},j_{2})$ for every $(i_{1},j_{1}),(i_{2},j_{2})\in\mathcal{E}(\mathcal{X})$, and that $\mathbb{E}\left[Y_{ij}^{[m]}\right]=d_{ij}^{\mathcal{X}}$ for every $m=1,\dots,M$. Under these assumptions, the likelihood can be expressed as, $f_{\mathscr{Y}}(\mathscr{y}\mid d_{\mathcal{E}(\mathcal{X})})=\prod_{(i,j)\in\mathcal{E}(\mathcal{X})}f_{\mathscr{Y}_{ij}}\left(\mathscr{y}_{ij}\bigm{|}d_{ij}^{\mathcal{X}}\right)$. Additionally, suppose that all measurements associated with edge $(i,j)$ are independent and identically distributed (i.i.d.). The likelihood can be further simplified,

$$f_{\mathscr{Y}}\left(\mathscr{y}\mid d_{\mathcal{E}(\mathcal{X})}\right)=\prod_{(i,j)\in\mathcal{E}(\mathcal{X})}\prod_{m=1}^{M}f_{Y_{ij}^{[m]}}\left(y_{ij}^{[m]}\Bigm{|}d_{ij}^{\mathcal{X}}\right),$$

(1)

where $f_{Y_{ij}^{[m]}}\left(y_{ij}^{[m]}\Bigm{|}d_{ij}^{\mathcal{X}}\right)$ is the true marginal pdf of $Y_{ij}^{[m]}$. In the rest of the paper, we will use $f_{Y_{ij}^{[m]}}\left(y_{ij}^{[m]}\Bigm{|}d_{ij}^{\mathcal{X}},\bm{\theta}_{ij}\right)$, where $\bm{\theta}_{ij}$ is a vector of shape parameters, to denote the marginal pdf for $Y_{ij}^{[m]}$. When a single shape parameter describes this marginal pdf, the scalar version $\theta_{ij}$ will be used instead.

Given a set of noisy edge length measurements $\mathscr{y}$, an ML estimate of $\mathcal{X}$ can be obtained by maximizing (1). This case will be referred to as matched DGP estimation. In practice, we often do not know what parametric form the marginal pdf of $Y_{ij}^{[m]}$ takes. In this case, we would modify  ​(1) by replacing the true marginal pdf of $Y_{ij}^{[m]}$ with an assumed version, $f_{Y_{ij}^{[m]}}^{*}\left(y_{ij}^{[m]}\Bigm{|}d_{ij}^{\mathcal{X}}\right)$. If $f_{Y_{ij}^{[m]}}^{*}\left(y_{ij}^{[m]}\Bigm{|}d_{ij}^{\mathcal{X}}\right)\neq f_{Y_{ij}^{[m]}}\left(y_{ij}^{[m]}\Bigm{|}d_{ij}^{\mathcal{X}}\right)$, then we refer to the modified version of the likelihood as the mismatched likelihood. Moreover, an estimate of $\mathcal{X}$ that is obtained by maximizing the mismatched likelihood will not necessarily be an ML estimate. We refer to this case as mismatched DGP estimation. We are interested in the difference in quality of estimates of $\mathcal{X}$ in the matched and mismatched cases.

Limitations of least-squares. Suppose that $Y_{ij}^{[m]}$ follows a Gaussian distribution with mean $d_{ij}^{\mathcal{X}}$ and variance $\theta_{ij}$, namely,

$$f_{Y_{ij}^{[m]}}\left(y_{ij}^{[m]}\Bigm{|}d_{ij}^{\mathcal{X}},\theta_{ij}\right)=\frac{1}{\sqrt{2\pi\theta_{ij}}}\exp\left(-\frac{(y_{ij}^{[m]}-d_{ij}^{\mathcal{X}})^{2}}{2\theta_{ij}}\right).$$

Then, maximizing the likelihood $f_{\mathscr{Y}}\left(\mathscr{y}\mid d_{\mathcal{E}(\mathcal{X})}\right)$ with respect to (w.r.t.) $\mathcal{X}$ is equivalent to minimizing the SSE w.r.t. $\mathcal{X}$. Several approaches to solve the DGP in the case of noisy measurements use variants of the SSE function [3, 8, 9]. A survey of prior art appears in Sec. II.

The aforementioned equivalence between minimizing the SSE function and maximizing the likelihood derived in (1) is only true when $Y_{ij}^{[m]}$ is Gaussian. In practice, depending on the measurement technique, it is unlikely that the noisy measurement $Y_{ij}^{[m]}$ is Gaussian distributed. Therefore, it is unlikely that minimizing the SSE function will yield an ML estimate of a structure.

Contributions. In this paper, we empirically show that, compared to mismatched DGP estimation, matched DGP estimation can offer an improvement of several dB. Moreover, our contribution is that, to the best of our knowledge, and in the context of the DGP, we are the first to focus on ML estimation beyond Gaussian-distributed measurements. The rest of the paper is structured as follows: we review the prior art in Sec. II. The different marginal pdfs of $Y_{ij}^{[m]}$ that will be considered, and how estimates of $\mathcal{X}$ will be evaluated, are discussed in Sec. III. Experimental results are presented in Sec. IV. Finally, we summarize our results and describe future work in Sec. V​.

One of the earliest discussions of cost functions related to the solution of the DGP was due to Crippen and Havel [3]. In their work, they used a variant of the SSE function, which we refer to as SSE-CH, to estimate molecular structures with physical constraints. Although others have empirically shown that SSE-CH yields acceptable results [3], the resulting estimate of a structure need not be an ML estimate.

Subsequent papers implement modifications to the SSE-CH cost function. For example, Moré and Wu [9] adopted a continuation approach: the SSE-CH cost function is repeatedly smoothed via a convolution with a Gaussian kernel. After each successive smoothing stage, a minimizer of the smoothed SSE-CH is obtained. Eventually, the minimizers are traced back to the original function [9, p. 827]. Although this approach yielded promising results, the underlying cost function that is smoothed is still the SSE function. Moreover, there is no guarantee that the resulting estimate of $\mathcal{X}$ is an ML estimate.

Souza et al. [10] attempt to make SSE-CH differentiable by approximating certain parts using hyperbolic functions. Souza et al. [11] made further improvements by applying these approximations to the geometric buildup algorithm proposed by Luo and Wu [12]. However, the underlying cost function is still a variant of the SSE function, and so there is no guarantee of obtaining ML estimates. For a detailed review of approaches to solve the DGP, see [7]. More recently, Gopalkrishnan et al. [13] optimize a variation of SSE-CH to reconstruct DNA molecules from pairwise distance measurements. Giamou et al. [14] also optimize a variation of SSE-CH to determine the joint angles of a robot given the robot’s forward kinematics and desired poses. However, because variations of SSE-CH are used, it is unlikely that estimates obtained in these applications are ML estimates.

We will consider different forms of the marginal pdf of $Y_{ij}^{[m]}$, as included in (1), that result in a likelihood $f_{\mathscr{Y}}\left(\mathscr{y}\mid d_{\mathcal{E}(\mathcal{X})}\right)$ that, when maximized w.r.t. $\mathcal{X}$, yields ML estimates. It was previously shown that if $Y_{ij}^{[m]}$ follows a Gaussian distribution, then minimizing the SSE function is equivalent to maximizing the likelihood in (1). We consider the cases when $Y_{ij}^{[m]}$ follows either a Laplace or non-standardized Student’s $t$ (NSST) distribution.

If $Y_{ij}^{[m]}$ follows a Laplace distribution with mean $d_{ij}^{\mathcal{X}}$ and variance $2\theta_{ij}^{2}$, then

$$f_{Y_{ij}^{[m]}}\left(y_{ij}^{[m]}\Bigm{|}d_{ij}^{\mathcal{X}},\theta_{ij}\right)=\frac{1}{2\theta_{ij}}e^{-\left\lvert y_{ij}^{[m]}-d_{ij}^{\mathcal{X}}\right\rvert\big{/}\theta_{ij}}.$$

(2)

If $Y_{ij}^{[m]}$ follows an NSST distribution with mean $d_{ij}^{\mathcal{X}}$ and variance $b_{ij}^{2}\nu_{ij}/(\nu_{ij}-2)$, assuming $\nu_{ij}>2$, then an expression for the marginal pdf of $Y_{ij}^{[m]}$ can be derived as follows. Let $W$ follow a Student’s $t$ distribution with parameter $\nu_{ij}$, which is a positive integer. The pdf of $W$ is

$$f_{W}(w\mid\nu_{ij})=\frac{\Gamma\left((\nu_{ij}+1)/2\right)}{\Gamma\left(\nu_{ij}/2\right)\sqrt{\nu_{ij}\pi}}\left(1+(w^{2}/\nu_{ij})\right)^{-(\nu_{ij}+1)/2},$$

where $\Gamma(z)=\int_{0}^{\infty}t^{z-1}e^{-t}\,dt$ for every $0<z\in\mathbb{R}$. If $Y_{ij}^{[m]}=b_{ij}W+d_{ij}^{\mathcal{X}}$, then

$$f_{Y_{ij}^{[m]}}\left(y_{ij}^{[m]}\Bigm{|}d_{ij}^{\mathcal{X}},\bm{\theta}_{ij}\right)=\frac{1}{b_{ij}}f_{W}\left(\frac{y_{ij}^{[m]}-d_{ij}^{\mathcal{X}}}{b_{ij}}\Biggm{|}\nu_{ij}\right),$$

(3)

where $\bm{\theta}_{ij}=(\nu_{ij},b_{ij})$.

Evaluation. Given an estimate $\widehat{\mathcal{X}}$ of $\mathcal{X}$, if we rotate, reflect, and/or translate the points in $\widehat{\mathcal{X}}$ by an equal amount, then $\widehat{\mathcal{X}}$ will still have the same edge lengths as $\mathcal{X}$. Therefore, a unique solution to the DGP cannot exist without additional constraints. Moreover, the proximity of $\widehat{\mathcal{X}}$ to $\mathcal{X}$ will need to be evaluated up to rotation, translation, and reflection. Let $\mathbf{X}$ and $\widehat{\mathbf{X}}$ be $K\times N$ matrices, where the $n$th columns of $\mathbf{X}$ and $\widehat{\mathbf{X}}$ correspond to the $n$th points in $\mathcal{X}$ and $\widehat{\mathcal{X}}$, respectively. We first translate the column vectors in $\mathbf{X}$ and $\widehat{\mathbf{X}}$ so that they are centered at the origin by computing $\mathbf{X}_{c}=\mathbf{X}-\mathbf{X}\mathbf{1}\big{/}N$ and $\widehat{\mathbf{X}}_{c}=\widehat{\mathbf{X}}-\widehat{\mathbf{X}}\mathbf{1}\big{/}N$, where $\mathbf{1}$ is an $N\times N$ matrix of ones.

Next, how close $\widehat{\mathcal{X}}$ is to $\mathcal{X}$ is equivalent to the solution of the orthogonal Procrustes problem (OPP) [15],

	$\displaystyle\min_{\mathbf{R}}$	$\displaystyle\left\|\mathbf{R}\widehat{\mathbf{X}}_{c}-\mathbf{X}_{c}\right\|_{F},$		(4)
	s.t.	$\displaystyle\mathbf{R}^{T}\mathbf{R}=\mathbf{I},$

where $\left\|\cdot\right\|_{F}$ is the Frobenius norm, which for a $K\times N$ matrix, $\mathbf{A}=(a_{ij})$, is defined as, $\left\|\mathbf{A}\right\|_{F}=\left(\sum_{k=1}^{K}\sum_{n=1}^{N}|a_{ij}|^{2}\right)^{1/2}$. The objective of the OPP is to determine the orthonormal matrix $\mathbf{R}$ that most closely maps $\widehat{\mathbf{X}}_{c}$ to $\mathbf{X}_{c}$. We will refer to the value obtained by solving the OPP for an estimate $\widehat{\mathcal{X}}$, divided by the number of points $N$, as the OPP loss associated with $\widehat{\mathcal{X}}$. As our objective is to show that an ML estimate $\widehat{\mathcal{X}}_{\text{ML}}$ of $\mathcal{X}$ more closely approximates $\mathcal{X}$ than a non-ML estimate $\widehat{\mathcal{X}}$, we will compare the OPP losses of $\widehat{\mathcal{X}}_{\text{ML}}$ and $\widehat{\mathcal{X}}$.

We begin our discussion with numerical results on triangles, which are structures consisting of 3 points (see Fig. 2). Later, we demonstrate analogous numerical results for structures with 10 points (see Fig. 5). In our experiments, we let $Y_{ij}^{[m]}$ follow either a Laplace (2) or NSST (3) distribution. If $Y_{ij}^{[m]}$ follows a Laplace (respectively, NSST) distribution, then the marginal pdf of $Y_{ij}^{[m]}$ in the likelihood (1) is varied between the Laplace (respectively, NSST) pdf in the matched case and the Gaussian pdf in the mismatched case.

The variance of the noisy measurements was changed as a function of different signal-to-noise ratio (SNR) values. If the distribution from which the i.i.d. noisy measurements associated with edge $(i,j)$ were obtained has variance $\sigma_{ij}^{2}$, and if we let $\sigma_{\mathcal{X}}$ be the average edge length of a structure $\mathcal{X}$, then $\text{SNR}=10\log_{10}\left(\sigma_{\mathcal{X}}^{2}/\sigma_{ij}^{2}\right)$, where SNR is quantified in decibel (dB) units. We leave the case where $\sigma_{ij}^{2}$ depends on edge $(i,j)$ for future work. The SNR was varied from –20 dB to 20 dB in steps of 5 dB. For each SNR value, we sample a set of noisy measurements, maximize the likelihood¹¹1An implementation of the L-BFGS algorithm [16] that is available in PyTorch [17] was used to maximize the likelihood. in (1) for the Laplace (or NSST) cases, and then report the OPP loss associated with an estimated structure $\widehat{\mathcal{X}}$ for each case.

Triangles. We first considered 8 triangles in 2D. Two of these triangles, together with example estimates and their associated OPP losses, are shown in Fig. 2. For each SNR value, and for each triangle, we sampled $M=10$ noisy measurements per edge 100 times. For each sample, we computed an estimate $\widehat{\mathcal{X}}$ using different likelihoods and then computed the OPP loss associated with $\widehat{\mathcal{X}}$. The distributions of these OPP losses are shown in Fig. 3(a) and Fig. 3(b). We also computed the pairwise difference in OPP losses resulting from mismatched and matched estimation. The distributions of these differences are shown in Fig. 3(c).

We observe that, in general, the distribution of OPP losses associated with the likelihood that is matched to the noise distribution is shifted lower than the distribution associated with the likelihood that is not matched to the noise distribution. Moreover, in several instances, the increase in median OPP loss due to a 3–4 dB decrease in SNR can approximately be compensated for by choosing to use the likelihood that is matched to the noise distribution.

We then repeated this experiment for different values of $M$, the number of noisy measurements per edge. The median OPP losses for $M=10,25,50,100$ are shown in Fig. 4. We observe that, in general, maximizing the likelihood that is matched to the noise distribution offers an improvement over mismatched DGP estimation comparable to the improvement obtained by doubling the number of noisy measurements per edge.

10 points. To further support our results, we repeated our experiments for 30 10-point structures in 2D. One of these structures is shown in Fig. 5. The distributions of OPP losses are shown in Fig. 6(a) and Fig. 6(b), and the distributions of pairwise differences in OPP losses for the mismatched and matched cases are shown in Fig. 6(c). Similar to the case of triangles, we observe that, in general, the distribution of OPP losses associated with the likelihood that is matched to the noise distribution is lower than the distribution associated with the likelihood that is not matched to the noise distribution.

In this paper, we have explored the impact of matching the likelihood function to the noise distribution on estimation accuracy in the context of the distance geometry problem (DGP). We have empirically found that matching the likelihood function to the noise distribution can offer a 3–4 dB compensation for the loss in estimation accuracy due to noisier measurements and smaller sample sizes. Moreover, matched estimation can be helpful in these scenarios, which often arise in practice.

There are several directions for future work that can be explored. For example, in practice, measurements may not be available for some edges, so the case where the numbers of edge measurements are unequal may also be considered. We may also explore the application of our approach to real-world measurements. Finally, a theoretical evaluation of the mismatch between the likelihood function and the noise distribution would be an insightful endeavor.

We thank Vaibhav Choudhary and Nikhil Gopalkrishnan for informative discussions.