# How can we objectively measure the similarity between two scatter plots whose coordinates are known?

I have two scatter plots $$S_1$$ showing the points $$(x_i, y_i)$$ and $$S_2$$ showing the points $$(u_i, v_i)$$ and both scatter plots are drawn in the same scale. Let us call $$S_1$$ as the base plot and $$S_2$$ as a estimated plot of the base. Thus for any $$i$$, the point $$(u_i, v_i)$$ is an estimate of the point $$(x_i, y_i)$$.

I want to objectively measure how closely the two scatter plots resemble each other such that :

• If $$S_2$$ is a identical to $$S_1$$ then the resemblance must be 100%
• If $$S_2$$ is an exact rotation $$S_1$$ around any point in the XY plane then also the resemblance must be 100%
• For all other cases the resemblance must between 0-100% depending how similar they are

A brute force approach is to rotate one plot in increments of 1 degree and use a convolution neural network based image recognition to measure the similarity between the two images or use other image recognition techniques such as SIFT.

However my use-case much simpler as I only have the coordinates of the points of the two plots so I am looking for a statistical / geometrical solution using only the knowledge of the $$(x_i, y_i)$$ and $$(u_i, v_i)$$ .

Question: How can we construct a metric that will measure the resemblance between two scatter plots?

• You've indexed the points in your two diagrams. Is that important? Once you apply a rotation, are you expecting $\left(x_1,y_1\right)$ to map to $\left(u_1,v_1\right)$, $\left(x_2,y_2\right)$ to $\left(u_2,v_2\right)$ and so on, or could the points be reordered? Commented Aug 14 at 8:39
• Do you want rotation around Origin ? Do you want rotation by arbitrary angle ? Do you additionally want mirror image ?
– Prem
Commented Aug 14 at 8:40
• @ChrisLewis Yes you are right. We can say that for each $i$, $(u_i,v_i)$ is estimate of the point $(x_i,y_i)$. I have edited the post to clarify this. Commented Aug 14 at 8:49
• One term for this is Procrustes analysis Commented Aug 14 at 22:41
• Check this question and its answer. Commented Aug 18 at 1:06

## 5 Answers

Let $$a_i=\begin{pmatrix} x_i \\ y_i\end{pmatrix}$$ and $$b_i=\begin{pmatrix} u_i \\ v_i\end{pmatrix}$$.

With a perfect match, we would have $$b_i = R(\theta)a_i+t$$

for every $$i$$, for some rotation matrix $$R(\theta)=\begin{pmatrix}\cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{pmatrix}$$ and translation vector $$t=\begin{pmatrix} t_1 \\ t_2 \end{pmatrix}$$.

It makes sense therefore to look at minimising

$$F(\theta,t)=\sum_i \left|R(\theta)a_i+t-b_i\right|^2$$

(a least-squares regression). Setting the partial derivative (gradient) with respect to $$t$$ to zero,

$$\sum_i R(\theta)a_i+t-b_i=0$$

ie

$$t=\frac1{n}\sum_i \left(b_i-R(\theta)a_i\right)$$

Likewise, setting the partial derivative with respect to $$\theta$$ to zero,

$$\sum_i R'(\theta)a_i \cdot \left(R(\theta)a_i+t-b_i\right)=0$$

Note that $$R'(\theta)=\begin{pmatrix}-\sin \theta & -\cos \theta \\ \cos \theta & -\sin \theta \end{pmatrix}$$ corresponds to a rotation through $$\theta+\frac{\pi}{2}$$, so $$R'(\theta)a_i$$ is always perpendicular to $$R(\theta)a_i$$. This reduces the last expression to

$$\sum_i R'(\theta)a_i \cdot \left(t-b_i\right)=0$$

Substituting for $$t$$ (and changing its summation index for clarity), and multiplying up by $$n$$,

$$\sum_i R'(\theta)a_i \cdot \left(\sum_j \left(b_j-R(\theta)a_j\right)-nb_i\right)=0\tag{*}$$

Now, I'm pretty sure there's a good linear algebra way to proceed here, but I can't find it, so my best suggestion at this point is to note that, although $$(*)$$ looks horrific, it's really just an equation in terms in $$\cos\theta$$, $$\sin\theta$$, $$\cos2\theta$$ and $$\sin2\theta$$. I don't think it will be solvable analytically, so it's probably simplest to calculate the coefficients based on the data at this stage and solve numerically. Once you have $$\theta$$, you can work out $$t$$ and then $$F$$ (which is the score you're after).

If anyone else has any ideas how to proceed from $$(*)$$ please let me know!

• There is a solution in terms of the singular value decomposition, see here for example. Commented Aug 14 at 22:44
• This has some good thinking , though I think it is unnecessary to calculate & minimize the expected rotation angle & the translation. OP just wants to know whether the Points have been rotated or not , the optimum rotation angle is not wanted.
– Prem
Commented Aug 16 at 15:43
• @Prem Thanks. The actual calculation seemed necessary in order to find the score (or a version of it) that OP was after. However, that was only the case because I couldn't get an analytical approach to work. I haven't quite got my head around the singular value decomposition yet; it may well be that that can yield a nice form for the score directly. Commented Aug 17 at 8:06
• Note that $R'(\theta)^TR(\theta)=R(-\frac\pi2)$, so that terms of the double angle should cancel. It remains something like $c=trace(R'(\theta)^TM)$. In the ideal case $c$ is a minimum of the right side, so that only one solution exists per period. Commented Aug 18 at 7:22

We want to know whether 2 "Polygons" $$S_1$$ & $$S_2$$ (which have the Vertex lists given via Scatter Plots) are same or not. That similarity is based on whether $$S_1$$ can be moved , shrunk or rotated to get $$S_2$$ , to generate a number between $$0 \%$$ & $$\pm 100 \%$$
These changes will not change the internal angles of the Polygons.

## Basic Algorithm :

Here is a very simple algorithm suiting that :
Let $$A_i$$ be the Points in $$S_1$$ & $$B_i$$ be the Points in $$S_2$$
Let $$A_1 A_2$$ be a reference vector which we consider to be rotated to $$B_1 B_2$$
We can calculate the angles $$\theta_i$$ between $$A_1 A_2$$ & $$A_1 A_i \text{ where } i=(2 \cdots n)$$
The Summation $$\Theta_1 = \Sigma \theta_i$$ (or the Average of that) of those angles gives some indication on how much the Polygon angles should add up to.
When $$S_1$$ is changed into $$S_2$$ , the changes should not change the internal angles of the Polygons.

[[ Image is indicating Angle-Pairs : the angle between $$A_1 A_2$$ & $$A_1 A_3$$ which should ideally be rotated to give the angle between $$B_1 B_2$$ & $$B_1 B_3$$ ; like-wise the angle between $$A_1 A_2$$ & $$A_1 A_6$$ which should ideally be comparable to the angle between $$B_1 B_2$$ & $$B_1 B_6$$ ]]

Hence we can compute the angles $$\theta_i$$ between $$B_1 B_2$$ & $$B_1 B_i \text{ where } i=(2 \cdots n)$$
We then can compute the Summation $$\Theta_2 = \Sigma \theta_i$$ (or the Average of that) of those angles.

Comparing $$\Theta_1$$ & $$\Theta_2$$ will tell us whether the Polygons are similar or not.
Let the larger magnitude value be $$100\%$$
Calculate the other value as a Percentage of that (including the Sign) to get the "Similarity Measure" we want.

## Some tweaking :

Comparing the total or average will make calculations very fast in linear time.
We can additionally calculate all Pairs of vectors $$A_1 A_n$$ & $$A_1 A_m$$ to get a larger total & then make the final calculation.
[[ That will make it quadratic time , though "Accuracy" is improved , in the sense that a few wrong angles might contribute less to the total ]]

Basic Algorithm gives the total angle to compare with.
Instead , we can generate the list of angles & then take the variance of those angles. Here $$0$$ variance indicates $$100\%$$ match.
[[ That will take more computation , though "Accuracy" might be more , in the sense that "Positive" rotations might not get canceled by "Negative" rotations ]]

There are more tweaks we can consider ( eg calculating the angles from the Centroid , calculating with mirror image ) , though listing those will make this Post unnecessarily long.

If $$a_i = [x_i, y_i]^T, b_i = [u_i, v_i]^T, i = 1, 2, \dots, N$$, then the error function we're trying to minimize is

$$E = \displaystyle \sum_{i=1}^N ( b_i - (t + s R a_i) )^T (b_i - (t+ s R a_i))$$

where $$s \gt 0$$ is a scale factor, and $$R$$ is a rotation matrix, and $$t$$ is a displacement vector.

In this setup, the sequence $$\{b_i\}$$ is allowed to be a scaled/rotated version of the sequence $$\{a_i\}$$ while still maintaining perfect similarity.

We want to find $$\theta$$ which determined the rotation matrix $$R$$, and $$s$$ that will result in minimum $$E$$.

Taking the partial derivative of $$E$$ with respect to $$t$$ gives

$$\nabla_t E = \displaystyle -2 \sum_{i=1}^N (b_i - (t + s R a_i)) = \mathbf{0}$$

Hence,

$$N t = \displaystyle \sum_{i=1}^N (b_i - s R a_i )$$

i.e.

$$t = \dfrac{1}{N} \displaystyle \sum_{i=1}^N (b_i - s R a_i )$$

Taking the partial derivative of $$E$$ with respect to $$\theta$$ gives

$$\dfrac{\partial E}{\partial \theta} = \displaystyle 2 \sum_{i=1}^N (b_i - (t + s R a_i))^T ( - \dfrac{\partial t }{\partial \theta} - s R' a_i ) = 0$$

This expands to

$$\displaystyle 2 \sum_{i=1}^N ( - b_i^T \dfrac{\partial t }{\partial \theta} - s b_i^T R' a_i + t^T \dfrac{\partial t }{\partial \theta} + s t^T R' a_i + s a_i^T R^T \dfrac{\partial t}{\partial \theta} + s^2 a_i^T R^T R' a_i) = 0$$

From the expression for $$t$$ above, we have

$$\dfrac{\partial t}{\partial \theta} = \dfrac{1}{N} \displaystyle \sum_{i=1}^N (-s R' a_i)$$

Therefore,

$$\displaystyle \sum_{i=1}^N b_i^T \dfrac{\partial t }{\partial \theta} = \dfrac{1}{N} \sum_{i = 1 }^N \sum_{j=1}^N b_i^T (-s R' a_j )$$

At this point we note that

$$R = \begin{bmatrix} \cos \theta && - \sin \theta \\ \sin \theta && \cos \theta \end{bmatrix} = \cos \theta \ I + \sin \theta \ Q$$

Where $$I$$ is identity matrix and $$Q = \begin{bmatrix} 0 && -1 \\ 1 && 0 \end{bmatrix}$$

Therefore

$$R' = \dfrac{\partial R}{\partial \theta} = - \sin \theta \ I + \cos \theta \ Q$$

And we also note that $$R^T R' = Q$$, and that $$Q$$ is skew-symmetric, so for any vector $$w$$, we have $$w^T Q w = 0$$.

Therefore,

$$\displaystyle \sum_{i=1}^N b_i^T \dfrac{\partial t }{\partial \theta} = \dfrac{-s}{N} \sum_{i = 1 }^N \sum_{j=1}^N b_i^T (- \sin \theta \ I + \cos \theta \ Q ) a_j$$

And this evaluates to

$$\displaystyle \sum_{i=1}^N b_i^T \dfrac{\partial t }{\partial \theta} = \dfrac{-s}{N} ( - \sin \theta \ S_b^T S_a + \cos \theta \ S_b^T Q S_a )$$

where $$S_a = \displaystyle \sum_{i=1}^N a_i$$ and $$S_b = \displaystyle \sum_{i=1}^N b_i$$

Moving on to the second term,

$$\displaystyle \sum_{i=1}^N b_i^T R' a_i = \sum_{i=1}^N b_i^T (- \sin \theta \ I + \cos \theta \ Q ) a_i = - \sin \theta S_1 + \cos \theta S_2$$

where $$S_1 = \displaystyle \sum_{i=1}^N b_i^T a_i , S_2 = \displaystyle \sum_{i=1}^N b_i^T Q a_i$$

Next, we have

$$\displaystyle \sum_{i=1}^N t^T \dfrac{\partial t }{\partial \theta} = \dfrac{1}{N^2} \sum_{i=1}^N \sum_{j=1}^N \sum_{k=1}^N (b_j - s R a_j)^T (- s R' a_k )$$

And this evaluates to

$$\displaystyle \sum_{i=1}^N t^T \dfrac{\partial t }{\partial \theta} = \dfrac{s}{N} ( \sin \theta \ S_b^T S_a - \cos \theta \ S_b^T Q S_a )$$

The next term is

$$\displaystyle \sum_{i=1}^N t^T R' a_i = \dfrac{1}{N} ( -\sin \theta S_b^T S_a + \cos \theta S_b^T Q S_a )$$

The next term is

$$\displaystyle \sum_{i=1}^N a_i^T R^T \dfrac{\partial t}{\partial \theta} = \dfrac{1}{N} \sum_{i=1}^N \sum_{j=1}^N - s a_i^T R^T R' a_j = 0$$

Also, the final terms is

$$\displaystyle \sum_{i=1}^N a_i^T R^T R' a_i = 0$$

This is because $$R^T R' = Q$$ which is a skew-symmetric matrix.

Now our equation becomes

$$\dfrac{s}{N} ( - \sin \theta \ S_b^T S_a + \cos \theta \ S_b^T Q S_a ) -s(- \sin \theta S_1 + \cos \theta S_2) + \dfrac{s}{N} ( \sin \theta \ S_b^T S_a - \cos \theta \ S_b^T Q S_a ) + \dfrac{s}{N} ( -\sin \theta S_b^T S_a + \cos \theta S_b^T Q S_a ) = 0$$

Dividing through by $$s$$, this equation becomes

$$\cos(\theta) ( \dfrac{1}{N} S_b^T Q S_a - S_2 ) + \sin(\theta) ( S_1 - \dfrac{1}{N} S_b^T S_a) = 0$$

Therefore, the optimal angle is given by

$$\boxed{\theta = \tan^{-1} \left( \dfrac{ N \ S_2 - S_b^T Q S_a }{ N \ S_1 - S_b^T S_a } \right)}$$

Similarly, we can find the partial derivative of $$E$$ with respect to $$s$$, and equate it to zero, and this gives

$$\dfrac{\partial E}{\partial s } = 2 \displaystyle \sum_{i=1}^N (b_i - (t + s R a_i))^T (- (\dfrac{\partial t}{\partial s } + R a_i ) ) = 0$$

We note that

$$\dfrac{\partial t}{\partial s} = \dfrac{1}{N} \displaystyle \sum_{i=1}^N (- R a_i) = - \dfrac{1}{N} R S_a$$

Expanding $$\dfrac{\partial E}{\partial s}$$ we get

$$\dfrac{\partial E}{\partial s} = 2 \displaystyle \sum_{i=1}^N (\dfrac{1}{N} b_i^T R S_a - b_i^T R a_i - \dfrac{1}{N} t^T R S_a + t^T R a_i - \frac{s}{N} a_i^T R^T R S_a + s a_i^T R^T R a_i) = 0$$

Recall that $$t = \dfrac{1}{N} (S_b - s R S_a)$$. Hence, the summation evaluates to,

$$\dfrac{1}{N} S_b^T R S_a - (\cos \theta S_1 + \sin \theta S_2) - \dfrac{1}{N} (S_b - s R S_a)^T R S_a + \dfrac{1}{N} (S_b - s R S_a)^T (R S_a) - \dfrac{s}{N} S_a^T S_a + s S_3$$

where $$S_3 = \displaystyle \sum_{i=1}^N a_i^T a_i$$

This reduces to

$$- (\cos \theta S_1 + \sin \theta S_2) + \dfrac{1}{N} (S_b^T R S_a) - \dfrac{s}{N} S_a^T S_a + s S_3 = 0$$

Substituting for $$R = \cos \theta \ I + \sin \theta \ Q$$, and multiplying through by $$N$$, gives

$$- N (\cos \theta S_1 + \sin \theta S_2) + (S_b^T (\cos \theta \ I + \sin \theta Q) S_a) - s S_a^T S_a + s N \ S_3 = 0$$

This equation is linear in $$s$$. Hence,

$$\boxed{s = \dfrac{ \cos(\theta) ( N S_1 - S_b^T S_a ) + \sin(\theta) ( N S_2 - S_b^T Q S_a) }{ N S_3 - S_a^T S_a } }$$

The two expressions for $$\theta$$ and $$s$$ derived above result in the minimum $$E$$. To find the quality of similarity between the two point sequences $$\{a_i\}$$ and $$\{b_i\}$$, we just have to evalute the error function $$E$$ at this these values of $$\theta$$ and $$s$$.

As a quantitative measure, we can take the root of the average error per point $$E$$, i.e. compute

$$\Delta = \sqrt{ \dfrac{E}{N} }$$

This gives an indication of how far (on average) is $$b_i$$ away from $$(t + s R a_i)$$.

You could define the similarity measure $$r = 1- \min_{Q^T Q = I} \frac{\|QA-B\|_F^2}{\|B\|_F^2}$$ where $$A,B$$ contain your respective data. This is zero if and only there is an orthogonal matrix that transforms $$A$$ to $$B$$.

The reason I am suggesting this is because there is a known closed form solution for it, called orthogonal Procrustes, that is easy to compute. Formally, we minimize $$\|Q A - B\|_F^2$$ over $$Q$$ orthonormal. This is the same as maximizing $$Trace(Q AB^T) = Trace(Q U\Lambda V^T) = Trace(V^T Q U \Lambda^T)$$ where $$U\Lambda V^T = AB^T$$ is the SVD. Since $$V^T, Q, U^T$$ are orthogonal, so is their product, which we denote $$R$$. Now maximizing $$Trace(R \Lambda)$$ over orthogonal $$P$$ is done if $$R=I$$, and therefore we simply let $$Q=VU^T$$.

• You're ignoring the shift between $A$ and $B$ that could be present. Commented Aug 18 at 23:14
• You're also ignoring scaling of $B$ relative to $A$. Commented Aug 18 at 23:17
• @that'swhatitis True, I was not sure whether that is what is wanted here. We can always normalize the data in both scatterplots to take mean shifts and scaling into account Commented Aug 19 at 7:06

I don't think that rating a resemblance of 0% makes sense and can be obtained by a formula.

What you need is one minus a ratio of the distance residue after rotation fitting (see other comments about the procrustes problem), to some reference distance.

For instance, the average distance between the corresponding points, over the average of the distances between all pairs of base points. Or the average distances from every point to its nearest neighbor. You can as well use the squared distance, giving more emphasis to the large ones.