机器学习--准备数据与Numpy(八)--距离矩阵计算

距离矩阵的计算

距离矩阵参考资料

在讲距离矩阵之前，先复习一下什么是 欧式距离 :
在做分类时，常常需要估算两个样本间的相似性度量（SimilarityMeasurement），这时经常就用到两个样本间的“距离”(Distance)，采用什么样的方法计算距离是很讲究，甚至关系到分类的正确与否。​
经常使用的度量方法是欧式距离，
欧氏距离是最易于理解的一种距离计算方法，源自欧氏空间中两点间的距离公式。

(1)二维平面上两点a(x1,y1)与b(x2,y2)间的欧氏距离：


(2)三维空间两点a(x1,y1,z1)与b(x2,y2,z2)间的欧氏距离：

(3)两个n维向量a(x11,x12,…,x1n)与 b(x21,x22,…,x2n)间的欧氏距离：

 也可以用表示成向量运算的形式：

例如：
k-nearest neighbor这部分中，核心就是计算训练集与测试集之元素之间的欧氏距离。要求从训练集取5000个图像，测试集取了500个图像，计 算这5000个用于训练的图像与500个用于测试的图像之间的欧氏距离，其结果就是一个5000*500的距离矩阵。

给定

阶矩阵

，满足

。这里第

列向量是

维向量。求

矩阵，使得：



这里提供4种方法，需要使用到以下Python库：

import numpy as np
import numpy.linalg as la

• 方法1：标准方法使用两层循环计算Dij

def compute_squared_EDM_method(X):
# determin dimensions of data matrix
m,n = X.shape
# initialize squared EDM D
D = np.zeros([n, n])
# iterate over upper triangle of D
for i in range(n):
for j in range(i+1, n):
D[i,j] = la.norm(X[:, i] - X[:, j])**2
D[j,i] = D[i,j]    #*1
return D

由于是计算矩阵自身行向量之间的距离，所以结果是一个对称的三角矩阵。注意*1行代码处所做的优化。

• 方法2：利用矩阵内积dot计算Dij

在第一种方法中，我们使用了numpy的norm这个方法，这个方法从数学上讲，其计算公式是：

然后我们又将这个计算结果平方后赋给

因此，我们实际上是在计算：

上述运算可以使用点积（即矩阵内积）来计算：

D[i,j] = np.dot(X[:,i]-X[:,j],(X[:,i]-X[:,j]).T)

def compute_squared_EDM_method2(X):
# determin dimensions of data matrix
m,n = X.shape
# initialize squared EDM D
D = np.zeros([n, n])
# iterate over upper triangle of D
for i in range(n):
for j in range(i+1, n):
d = X[:,i] - X[:,j]
D[i,j] = np.dot(d, d)
D[j,i] = D[i,j]
return D

• 方法3：避免循环内点积运算，减少dot调用次数

注意在上面的方法中，dot运算被调用了

次，并且每次进行了

次乘积运算和

次加法运算。尽管numpy底层可能对点积运算做了优化，但这里还是存在可能进行进一步优化。请看下面的数学推导：


这里

属于格拉姆矩阵中的元素，可以通过在循环外计算矩阵，在循环内直接引用元素值即可，从而在循环内我们只需要做两次加（减）法运算：

格拉姆矩阵的求法很简单，只需要：

现在代码变为：

def compute_squared_EDM_method3(X):
# determin dimensions of data matrix
m,n = X.shape
# compute Gram matrix
G = np.dot(X.T, X)
# initialize squared EDM D
D = np.zeros([n, n])
# iterate over upper triangle of D
for i in range(n):
for j in range(i+1, n):
d = X[:,i] - X[:,j]
D[i,j] = G[i,i] - 2 * G[i,j] + G[j,j]
D[j,i] = D[i,j]
return D

• 方法4：利用重复操作替代外部循环

假设距离矩阵可以表示为

，与公式

进行对比，有：

这里H中第i行的每一个元素，取值都为

，也就是H的每一列，都对应着格拉姆矩阵的对角阵，因此，我们可以用下面的代码来计算H：

H = np.tile(np.diag(G), (n,1))

此外，由于

,所以最终距离矩阵可以计算为

现在，代码不再需要循环了：

def compute_squared_EDM_method4(x):
m,n = X.shape
G = np.dot(X.T, X)
H = np.tile(np.diag(G), (n,1))return H + H.T - 2*G

# -*- coding:utf-8 -*-
# Author:Justin Chan
# Python:3.6

import numpy as np
import numpy.random as np_randon

X = np_randon.randint(10,size=(4,3))
print(X)

print(X[:,1])#矩阵第一列向量
print(X[1,:])#矩阵第一行向量

[[5 8 2]
[7 4 8]
[6 7 1]
[2 2 0]]
[8 4 7 2]
[7 4 8]

import numpy as np
import numpy.linalg as la
import time

X = np.array([range(0, 500), range(500, 1000)])
m, n = X.shape

#方法1：标准方法使用两层循环计算Dij
#时间复杂度O(n*n)
t = time.time()
D = np.zeros([n, n])
for i in range(n):
for j in range(i + 1, n):
D[i, j] = la.norm(X[:, i] - X[:, j]) ** 2
D[j, i] = D[i, j]
print(time.time() - t)

#方法2：利用矩阵内积dot计算Dij
# 时间复杂度O(n*n)*O(m)
t = time.time()
D = np.zeros([n, n])
for i in range(n):
for j in range(i + 1, n):
d = X[:, i] - X[:, j]
D[i, j] = np.dot(d, d)
D[j, i] = D[i, j]
print(time.time() - t)

#方法3：避免循环内点积运算，减少dot调用次数
# 时间复杂度O(n*n)
t = time.time()
G = np.dot(X.T, X)
D = np.zeros([n, n])
for i in range(n):
for j in range(i + 1, n):
D[i, j] = G[i, i] - G[i, j] * 2 + G[j,j]
D[j, i] = D[i, j]
print(time.time() - t)

#方法4：利用重复操作替代外部循环
t = time.time()
G = np.dot(X.T, X)
#把G对角线元素拎出来，列不变，行复制n遍。
H = np.tile(np.diag(G), (n, 1))
D = H + H.T - G * 2
print(time.time() - t)

1.051072120666504
0.3479161262512207
0.15077900886535645
0.003982067108154297

其它习题：

1.将numpy数组保存到文件并读取。
2.用numpy解一个线性方程组。
3.自己实现一遍例题。