How to implement finding the coordinates of points from distance matrix in python based on gram-matrix?
I want to study a bus stop optimization problem. However, I'm now stuck in how to convert the distance matrix to the real coordinates of points.
I have browsed a lot resouce and known using the formula：M(i, j) = 0.5(D(1, j)^2 + D(i, 1)^2 - D(i, j)^2)* to solve the problem enter link description here. I'm not good at math and I just want to implement to it.
Firstly, I try to understand the principle of Mathematics and here is my solution. enter link description here.
Then, I want to implment the algorithm using python for the following example. Here is my matrix which represents different distance for each bus-station. I want to transfer it to the coordinates of points.
Here is my code for implmenting:
import csv
import numpy as np
import math
class csv_util():
def generate_coordinate_point(self):
'''transfer the distance matrix to the real coordinate points'''
sqrt_result = 2*math.sqrt(2)
matrix = np.array([[0,2,2,sqrt_result],[2,0,sqrt_result,2],[2,sqrt_result,0,2],[sqrt_result,2,2,0]])
gram_matrix = self.calculate_gram_matrix(matrix)
a, b = np.linalg.eig(gram_matrix)
#b = b.astype(np.int16)
a = a.astype(np.int16)
eigen_vector = format(b)
length = a.size
tmp_matrix = np.zeros(length * length)
random_point_matrix = tmp_matrix.reshape(length, length)
for item1 in range(length):
random_point_matrix[item1][item1] = a[item1]
print("the eigen-value is: " + format(random_point_matrix))
print("the eigen-vector is: " + eigen_vector)
new_matrix = (np.sqrt(random_point_matrix))*b
print("the coordinate points: "+format(new_matrix))
def calculate_gram_matrix(self,matrix):
'''get the gram matrix for transfer to the coordinate points'''
length = matrix[0].size
tmp_matrix = np.zeros(length*length)
gram_matrix = tmp_matrix.reshape(length,length)
for item1 in range(length):
for item2 in range(length):
gram_matrix[item1][item2] = (math.pow(matrix[0][item2],2)+math.pow(matrix[0][item1],2)-math.pow(matrix[item1][item2],2))/2
if gram_matrix[item1][item2]<0.1 and gram_matrix[item1][item2]>-0.1:
gram_matrix[item1][item2] = 0
return gram_matrix
However, the result of final matrix isn't correct. The result like this:
the eigen-value is: [[12. 0. 0. 0.]
[ 0. 0. 0. 0.]
[ 0. 0. 4. 0.]
[ 0. 0. 0. 0.]]
-------------
the eigen-vector is: [[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00]
[ 4.08248290e-01 -5.77350269e-01 -7.07106781e-01 0.00000000e+00]
[ 4.08248290e-01 -5.77350269e-01 7.07106781e-01 0.00000000e+00]
[ 8.16496581e-01 5.77350269e-01 1.57009246e-16 0.00000000e+00]]
-------------
the coordinate points: [[ 0. 0. 0. 0. ]
[ 0. -0. -0. 0. ]
[ 0. -0. 1.41421356 0. ]
[ 0. 0. 0. 0. ]]
The final points like this: [0,0],[-0.0,-0.0],[-0.0,1.414421],[0.0,0.0]. They cannot be satisfied with the distance matrix in this example. Please help me how to get the correct points. Thanks!
The construction of the Gram matrix of dot products associated to the distance matrix and its further factorization is generally a great method, which also allows you to infer the dimension of the coordinate realization of the distance matrix. However, if in your case the realization is planar (two dimensional), then I think it is (arguably) easier (and possibly faster) to just approach it a bit more geometrically (again, you should be certain that the distance matrix is for points in 2D):
import numpy as np
import math
def x_coord_of_point(D, j):
return ( D[0,j]**2 + D[0,1]**2 - D[1,j]**2 ) / ( 2*D[0,1] )
def coords_of_point(D, j):
x = x_coord_of_point(D, j)
return np.array([x, math.sqrt( D[0,j]**2 - x**2 )])
def calculate_positions(D):
(m, n) = D.shape
P = np.zeros( (n, 2) )
tr = ( min(min(D[2,0:2]), min(D[2,3:n])) / 2)**2
P[1,0] = D[0,1]
P[2,:] = coords_of_point(D, 2)
for j in range(3,n):
P[j,:] = coords_of_point(D, j)
if abs( np.dot(P[j,:] - P[2,:], P[j,:] - P[2,:]) - D[2,j]**2 ) > tr:
P[j,1] = - P[j,1]
return P
sqrt_result = 2*math.sqrt(2)
D = np.array([[0, 2, 2, sqrt_result],
[2, 0, sqrt_result, 2],
[2, sqrt_result, 0, 2],
[sqrt_result, 2, 2, 0]])
P = calculate_positions(D)
print(P)
You may want to add some checks and improvements to make sure that the vectors P[1,:] and P[2,:] are not aligned, which is equivalent to checking that
abs( P[1,0]*P[2,1] - P[1,1]*P[2,0] ) < 0.0001 (or some more appropriate threshold)
If they are, simply implement a while loop until you find the first vector P[j0, :] that is not aligned with P[1,0]. The role of this first vector P[j0,:] unaligned with the initial vector P[1,:] allows you to have a useful if clause in the function vector(D). I did not include it to avoid obscuring the idea of the code.