文章目录
1. 引言
距离矩阵是一个包含一组点两两距离的矩阵(即 二维数组)。因此给定
N
N
N个欧几里得空间中的点,其距离矩阵就是一个非负实数作为元素的
N
×
N
N \times N
N×N的对称矩阵。在机器学习中距离矩阵都计算非常常见(只要涉及距离计算,基本都需要计算距离矩阵),在本篇博客中就来记录一下如何使用Python
都科学计算包numpy
计算向量都距离矩阵。本篇博客讲解以行向量的欧氏距离为例讲解,但是同时给出了列向量的代码,距离矩阵的数学表达为:
给定 n × m n \times m n×m阶矩阵 X X X,满足 X = [ x 1 x 2 . . . x n ] X= \begin{bmatrix} x_1\\x_2\\...\\x_n\end{bmatrix} X=⎣⎢⎢⎡x1x2...xn⎦⎥⎥⎤,里第 i i i行向量是 m m m维向量**(注意这里是行向量)**。使得:
D i j = ∥ x i − x j ∥ D_{ij}=\left \|x_i-x_j \right \| Dij=∥xi−xj∥
2. 矩阵自身(行)向量之间的距离矩阵计算
这里提供4种方法,需要使用到以下Python库:
import numpy as np
import numpy.linalg as la #和线性代数相关的库
2.1 第一种方法:简单使用两重循环
# 行向量
#[[1,2,3],
# [4,5,6]]
# 得到
#[[ 0, 5.196]
# [ 5.196, 0 ]]
def compute_squared_EDM_method(X):
# 获得矩阵都行和列,因为是行向量,因此一共有n个向量
n,m = X.shape
# 因为有n个向量,距离矩阵是n x n
D = np.zeros([n, n])
# 迭代求解向量的距离
for i in range(n):
for j in range(i+1, n):
# la.norm()求向量都范数,默认是2范数
D[i,j] = la.norm(X[i, :] - X[j, :])
D[j,i] = D[i,j]
return D
# 列向量
#[[1,2,3],
# [4,5,6]]
# 得到
#[[0, 1.414, 2.828]
# [1.414, 0, 1.414]
# [2.828, 1.414, 0 ]]
def compute_squared_EDM_method(X):
# 获得矩阵都行和列,因为是列向量,因此一共有m个向量
n,m = X.shape
# 因为有m个向量,距离矩阵是m x m
D = np.zeros([m, m])
# 迭代求解向量的距离
for i in range(m):
for j in range(i+1, m):
# la.norm()求向量都范数,默认是2范数(注意这里是列向量)
D[i,j] = la.norm(X[:, i] - X[:, j])
D[j,i] = D[i,j] #*1
return D
由于是计算矩阵自身向量之间的距离,所以结果是一个对称的三角矩阵。注意*1行代码处所做的优化。在上述方法中我们使用了两层循环,因此代码虽不简洁,但十分易懂。
2.2 第二种方法:矩阵內积双重循环
在第一种方法中,我们使用了numpy
的norm
这个方法,这个方法从数学上讲,其计算公式是:
∥
x
i
−
x
j
∥
=
(
x
i
−
x
j
)
(
x
i
−
x
j
)
T
\left \|x_i-x_j \right \| =\sqrt{(x_i-x_j)(x_i-x_j)^T}
∥xi−xj∥=(xi−xj)(xi−xj)T
但是从另一方面来讲,我们可以先求点积运算,然后在进行求根运算:
D
i
j
=
(
x
i
−
x
j
)
(
x
i
−
x
j
)
T
D_{ij}=\sqrt{(x_i-x_j)(x_i-x_j)^T}
Dij=(xi−xj)(xi−xj)T
上述运算可以使用点积(即矩阵内积)来计算:
# 这里是列向量
D[i,j] = np.sqrt(np.dot(X[:,i]-X[:,j],(X[:,i]-X[:,j]).T))
现在代码变化为:
# 行向量
#[[1,2,3],
# [4,5,6]]
# 得到
#[[ 0, 5.196]
# [ 5.196, 0 ]]
def compute_squared_EDM_method2(X):
# 获得矩阵都行和列,因为是行向量,因此一共有n个向量
n,m = X.shape
# 因为有n个向量,距离矩阵是n x n
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.sqrt(np.dot(d, d))
D[j,i] = D[i,j]
return D
# 列向量
#[[1,2,3],
# [4,5,6]]
# 得到
#[[0, 1.414, 2.828]
# [1.414, 0, 1.414]
# [2.828, 1.414, 0 ]]
def compute_squared_EDM_method2(X):
# 获得矩阵都行和列,因为是列向量,因此一共有m个向量
n,m = X.shape
# 因为有m个向量,距离矩阵是m x m
D = np.zeros([m, m])
# 迭代求解向量的距离
for i in range(m):
for j in range(i+1, m):
# 因为是列向量,这里是列索引
d = X[:,i] - X[:,j]
# 向量內积运算,并求根运算
D[i,j] = np.sqrt(np.dot(d, d))
D[j,i] = D[i,j]
return D
2.3 第三种方法:避免循环内的点积运算
注意在上面的方法中,dot
运算被调用了
n
2
−
n
n^2-n
n2−n(针对列向量,如果是行向量就是
m
2
−
m
m^2-m
m2−m)次,并且每次进行了m
(针对列向量,如果是行向量就是
n
n
n)次乘积运算和n-1
(针对列向量,如果是行向量就是
m
2
−
m
m^2-m
m2−m)次加法运算。尽管numpy
底层可能对点积运算做了优化,但这里还是存在可能进行进一步优化。请看下面的数学推导(行向量):
D i j = ( x i − x j ) ( x i − x j ) T = x i x i T − 2 x i x j T + x j x j T D_{ij} =\sqrt{(x_i-x_j)(x_i-x_j)^T}=\sqrt{x_ix_i^T-2x_ix_j^T+x_jx_j^T} Dij=(xi−xj)(xi−xj)T=xixiT−2xixjT+xjxjT
这里 x i x i T , x j x j T , x i x j T x_ix_i^T,x_jx_j^T,x_ix_j^T xixiT,xjxjT,xixjT属于格拉姆矩阵中的元素,可以通过在循环外计算矩阵,在循环内直接引用元素值即可,从而在循环内我们只需要做两次加(减)法运算:
D i j = G i i − 2 G i j + G j j D_{ij}=\sqrt{G_{ii}-2G_{ij}+G_{jj}} Dij=Gii−2Gij+Gjj
格拉姆矩阵的求法很简单,只需要:
# 行向量
G=np.dot(X, X.T)
# 列向量
G=np.dot(X.T, X)
现在代码变为:
# 行向量
# 行向量
#[[1,2,3],
# [4,5,6]]
# 得到
#[[ 0, 5.196]
# [ 5.196, 0 ]]
def compute_squared_EDM_method3(X):
# 获得矩阵的行和列,因为是行向量,因此一共有n个向量
n,m = X.shape
# 计算Gram 矩阵
G = np.dot(X, X.T)
# 初始化距离矩阵,因为有n个向量,距离矩阵是n x n
D = np.zeros([n, n])
# 迭代求解
for i in range(n):
for j in range(i+1, n):
D[i,j] = np.sqrt(G[i,i] - 2 * G[i,j] + G[j,j])
D[j,i] = D[i,j]
return D
# 列向量
#[[1,2,3],
# [4,5,6]]
# 得到
#[[0, 1.414, 2.828]
# [1.414, 0, 1.414]
# [2.828, 1.414, 0 ]]
def compute_squared_EDM_method3(X):
# 获得矩阵都行和列,因为是列向量,因此一共有m个向量
n,m = X.shape
# 计算Gram 矩阵
G = np.dot(X.T, X)
# 初始化距离矩阵, # 因为有m个向量,距离矩阵是m x m
D = np.zeros([m, m])
# 迭代求解
for i in range(m):
for j in range(i+1, m):
D[i,j] = np.sqrt(G[i,i] - 2 * G[i,j] + G[j,j])
D[j,i] = D[i,j]
return D
2.4 第四种方法:避免循环
假设距离矩阵可以表示为 D = H + K − 2 G D = \sqrt{H+K-2G} D=H+K−2G,与公式 D i j = G i i − 2 G i j + G j j D_{ij}=\sqrt{G_{ii}-2G_{ij}+G_{jj}} Dij=Gii−2Gij+Gjj进行对比,有:
H
i
j
=
G
i
i
,
K
i
j
=
G
j
j
H_{ij}=G_{ii}, K_{ij}=G_{jj}
Hij=Gii,Kij=Gjj
这里H
中第i
行的每一个元素,取值都为
G
i
i
G_{ii}
Gii,也就是H
的每一列,都对应着格拉姆矩阵的对角阵,因此,我们可以用下面的代码来计算H
(n是向量的个数,无论行向量还是列向量)(帮助理解:将G都对角线元素展开成一个行向量,那么
H
i
j
H_{ij}
Hij的每一个元素只和i有关,即相应G位置向上找,所以将行向量数值排列):
H
=
n
p
.
t
i
l
e
(
n
p
.
d
i
a
g
(
G
)
,
(
n
,
1
)
)
(将Y轴复制n倍,将X轴复制1倍)
H = np.tile(np.diag(G), (n,1)) \tag{将Y轴复制n倍,将X轴复制1倍}
H=np.tile(np.diag(G),(n,1))(将Y轴复制n倍,将X轴复制1倍)
此外,由于
K
=
H
T
K= H^T
K=HT,所以最终距离矩阵可以计算为
D
=
H
+
H
T
−
2
G
(开根代表对矩阵每一个元素求根)
D=\sqrt{H+H^T-2G} \tag{开根代表对矩阵每一个元素求根}
D=H+HT−2G(开根代表对矩阵每一个元素求根)
现在,代码不再需要循环了:
# 行向量
#[[1,2,3],
# [4,5,6]]
# 得到
#[[ 0, 5.196]
# [ 5.196, 0 ]]
def compute_squared_EDM_method4(X):
# 获得矩阵都行和列,因为是行向量,因此一共有n个向量
n,m = X.shape
# 计算Gram 矩阵
G = np.dot(X,X.T)
# 因为是行向量,n是向量个数,沿y轴复制n倍,x轴复制一倍
H = np.tile(np.diag(G), (n,1))
return np.sqrt(H + H.T - 2*G)
# 列向量
#[[1,2,3],
# [4,5,6]]
# 得到
#[[0, 1.414, 2.828]
# [1.414, 0, 1.414]
# [2.828, 1.414, 0 ]]
def compute_squared_EDM_method4(X):
# 获得矩阵都行和列,因为是列向量,因此一共有m个向量
n,m = X.shape
# 计算Gram 矩阵
G = np.dot(X.T, X)
# 因为是列向量,n是向量个数,沿y轴复制m倍,x轴复制一倍
H = np.tile(np.diag(G), (m,1))
return np.sqrt(H + H.T - 2*G)
2.5 第五种方法:利用scipy
求距离矩阵(推荐用法)
在scipy
中提供了一个工具,用于求距离矩阵,但是此工具算出来都结果不是一个矩阵,而是一个列表,此列表为距离矩阵的上三角的排列展开,如果想得到矩阵,需要转换,代码如下所示:
# 默认是针对行向量进行操作
# 向量矩阵为:
# [[1,2],
# [3,4],
# [5,6]
# [7,8]]
# 距离矩阵为:
# [[0, 2.828, 5.656, 8.485],
# [2.828, 0, 2.828, 5.656],
# [5.656, 2.828, 0, 2.828],
# [8.485, 5.656, 2.828, 0 ]]
# distA距离列表为(上三角矩阵展开成一个列表):
# [2.828, 5.656, 8.485, 2.828, 5.656, 2.828]
# distB距离矩阵为:
# [[0, 2.828, 5.656, 8.485],
# [2.828, 0, 2.828, 5.656],
# [5.656, 2.828, 0, 2.828],
# [8.485, 5.656, 2.828, 0 ]]
import numpy as np
from scipy.spatial.distance import pdist
from scipy.spatial.distance import squareform
A=np.array([[1,2],
[3,4],
[5,6],
[7,8]])
# A是一个向量矩阵:euclidean代表欧式距离
distA=pdist(A,metric='euclidean')
# 将distA数组变成一个矩阵
distB = squareform(distA)
3.两个矩阵之间的距离矩阵计算
3.1 第一种方法:使用numpy
计算
给定训练矩阵A为
5000
×
3072
5000\times 3072
5000×3072阶矩阵。这里5000
代表5000
幅带标签的图,3072
是其各像素在RGB
三个通道下的取值数。给定测试集矩阵B
为
500
×
3072
500\times 3072
500×3072阶矩阵。求矩阵B
的各行与矩阵A
的各行的距离(即两幅图的差异)矩阵,这个矩阵是一个
5000
×
500
5000\times 500
5000×500的矩阵。
更一般地,这个问题可以描述如下:
给定矩阵A为
m
×
k
m\times k
m×k阶矩阵,矩阵B为
n
×
k
n\times k
n×k阶矩阵,求矩阵B的任意行向量与矩阵A的任意行向量的距离矩阵
D
m
×
n
D_{m\times n}
Dm×n。这个矩阵的数学表达式为(
a
,
b
a, b
a,b均为行向量):
D
i
j
=
(
a
i
−
b
j
)
(
a
i
−
b
j
)
T
=
a
i
a
i
T
−
a
i
b
j
T
−
b
j
a
i
T
+
b
j
b
j
T
D_{ij}=\sqrt{(a_i-b_j)(a_i-b_j)^T}=\sqrt{a_ia^T_i-a_ib^T_j-b_ja^T_i+b_jb^T_j}
Dij=(ai−bj)(ai−bj)T=aiaiT−aibjT−bjaiT+bjbjT
为方便讨论,我们将上述各项分别记为
H
,
M
,
N
,
K
H, M, N, K
H,M,N,K,即:
D
m
×
n
=
H
m
×
m
−
M
m
×
n
−
N
n
×
m
+
K
n
×
n
(开根代表对矩阵每一个元素求根)
D_{m\times n}=\sqrt{H_{m\times m}-M_{m\times n}-N_{n\times m}+K_{n\times n}} \tag{开根代表对矩阵每一个元素求根}
Dm×n=Hm×m−Mm×n−Nn×m+Kn×n(开根代表对矩阵每一个元素求根)
显然上述公式是无法进行运算的,因为除了M与D外,其它矩阵的秩各不相同。所以我们要回到前一个数学表达式上。
- H H H对 D D D的贡献是对于 D D D的每一行,都加上 a i a i T a_ia^T_i aiaiT
- K K K对 D D D的贡献是对于 D D D的每一列,都加上 b j b j T b_jb^T_j bjbjT
- M M M和 N N N互为转置矩阵。即对 D i j D_{ij} Dij,要减去矩阵 N j , i N_{j,i} Nj,i元素,而这个元素就是 M i j M_{ij} Mij
- 由前三个总结 D i j D_{ij} Dij为 a i a i T − a i b j T − b j a i T + b j b j T {a_ia^T_i-a_ib^T_j-b_ja^T_i+b_jb^T_j} aiaiT−aibjT−bjaiT+bjbjT求根
代码如下(行向量):
# 行向量:A (3行2列)
#[[1,2],
# [3,4],
# [5,6]]
# 行向量:B (2行2列)
#[[1,2],
# [3,4]]
#得到矩阵C(3行2列),由A->B的距离,Cij代表A中都第i个向量到B中第j向量都距离
#[[0, 2.828],
# [2.828 , 0 ],
# [5.656 , 2.828]]
def compute_distances_no_loops(A, B):
#A 有m个向量
m = np.shape(A)[0]
#B 有n个向量
n = np.shape(B)[0]
# 求得矩阵M为 m*n维(针对行向量)
M = np.dot(A, B.T)
# 对于H,我们只需要A . A^T的对角线元素
# np.square(A)是A中都每一个元素都求平方
# np.square(A).sum(axis=1) 是将每一行都元素都求和,axis是按行求和(原因是行向量)
# np.matrix() 是将一个列表转为矩阵,该矩阵为一行多列
# 求矩阵都转置,为了变成一列多行
# np.tile是复制,沿Y轴复制1倍(相当于没有复制),再沿X轴复制n倍
H = np.tile(np.matrix(np.square(A).sum(axis=1)).T,(1,n))
# 对于H,我们只需要B . B^T的对角线元素
# np.square(B)是B中都每一个元素都求平方
# np.square(B).sum(axis=1) 是将每一行都元素都求和,axis是按行求和(原因是行向量)
# np.matrix() 是将一个列表转为矩阵,该矩阵为一行多列
# np.tile是复制,沿Y轴复制m倍(相当于没有复制),再沿X轴复制1倍
K = np.tile(np.matrix(np.square(B).sum(axis=1)),(m,1))
# H对M在y轴方向上传播,即H加和到M上的第一行,K对M在x轴方向上传播,即K加和到M上的每一列
return np.sqrt(-2 * M + H + K)
代码如下(列向量):
# 行向量:A (2行3列).3个向量
#[[1,2,3],
# [4,5,6]]
# 行向量:B (2行2列),2个向量
#[[1,2],
# [3,4]]
#得到矩阵C(3行2列),由A->B都距离 Cij代表A中都第i个向量到B中第j向量都距离
#[[1 , 1 ],
# [2.236, 1 ],
# [3.605, 2.236]]
def compute_distances_no_loops(A, B):
#A 有m个向量(针对列向量)
m = np.shape(A)[1]
#B 有n个向量(针对列向量)
n = np.shape(B)[1]
# 求得矩阵M为 m*n维
# 求得矩阵M为 m*n维
M = np.dot(A.T, B)
# 对于H,我们只需要A . A^T的对角线元素,下面的方法高效求解(只计算对角线元素)
#沿Y轴复制1倍(相当于没有复制),再沿X轴复制n倍
H = np.tile(np.matrix(np.square(A).sum(axis=0)).T,(1,n))
# 结果K为n维行向量.要将其元素运用到矩阵M的每一列,需要将其转置为行向量
K = np.tile(np.matrix(np.square(B).sum(axis=0)),(m,1))
# H对M在y轴方向上传播,即H加和到M上的第一行,K对M在x轴方向上传播,即K加和到M上的每一列
return np.sqrt(-2 * M + H + K)
3.2 第二种方法:利用scipy
求距离矩阵(推荐用法)
在scipy
中提供了一个工具,用于求两个向量集合距离矩阵,代码如下所示:
# 行向量:A (3行2列)
#[[1,2],
# [3,4],
# [5,6]]
# 行向量:B (2行2列)
#[[1,2],
# [3,4]]
#得到矩阵C(3行2列),由A->B的距离,Cij代表A中都第i个向量到B中第j向量都距离
#[[0, 2.828],
# [2.828 , 0 ],
# [5.656 , 2.828]]
import numpy as np
from scipy.spatial.distance import cdist
A=np.array([[1,2],
[3,4],
[5,6]])
B=np.array([[1,2],
[3,4]])
dist=cdist(A,B,metric='euclidean')
4. 相关数学知识:什么是格拉姆矩阵?
GRAM
中文名称为格拉姆矩阵,它是个有广泛应用的矩阵。
-
v
1
,
v
2
,
⋯
,
v
n
v_1,v_2,⋯,v_n
v1,v2,⋯,vn 是内积空间的一组行(列)向量,
Gram
矩阵定义为: G i j = ⟨ v i , v j ⟩ = v i v j T G_{ij}=⟨v_i,v_j⟩=v_iv_j^T Gij=⟨vi,vj⟩=vivjT,显然其是对称矩阵。 - 其实对于一个
X
N
⋅
d
X_{N⋅d}
XN⋅d(行向量
N
个样本,d
个属性)的样本矩阵而言, X ⋅ X T X⋅X^T X⋅XT 即为Gram
矩阵。如果是列向量 X T ⋅ X X^T⋅X XT⋅X 为Gram
矩阵(很重要,这里需要用到的) - 如果
v
1
,
v
2
,
⋯
,
v
n
v_1,v_2,⋯,v_n
v1,v2,⋯,vn分别是随机向量,则
Gram
矩阵是协方差矩阵; - 欧式空间中向量
v
1
,
v
2
,
⋯
,
v
n
v_1,v_2,⋯,v_n
v1,v2,⋯,vn的
Gram
矩阵一定是半正定矩阵,是正定矩阵的充要条件是 v 1 , v 2 , ⋯ , v n v_1,v_2,⋯,v_n v1,v2,⋯,vn线性无关。
具体形式为:
n
n
n维欧式空间中任意
k
(
k
≤
n
)
k(k≤n)
k(k≤n)个向量
α
1
,
α
2
,
⋯
,
α
k
α_1,α_2,⋯,α_k
α1,α2,⋯,αk的内积所组成的矩阵
Δ
(
α
1
,
α
2
,
⋯
,
α
k
)
=
(
(
α
1
,
α
1
)
(
α
1
,
α
2
)
.
.
.
(
α
1
,
α
k
)
(
α
2
,
α
1
)
(
α
2
,
α
2
)
.
.
.
(
α
2
,
α
k
)
.
.
.
.
.
.
.
.
.
.
.
.
(
α
k
,
α
1
)
(
α
k
,
α
2
)
.
.
.
(
α
k
,
α
k
)
)
\Delta\left(α_1,α_2,⋯,α_k\right)= \left(\begin{matrix} (α_1,α_1) & (α_1,α_2) & ... &(α_1,α_k) \\ (α_2,α_1) & (α_2,α_2) & ... & (α_2,α_k)\\ ... &...& ... &...\\ (α_k,α_1) & (α_k,α_2) & ... & (α_k,α_k) \end{matrix} \right)
Δ(α1,α2,⋯,αk)=⎝⎜⎜⎛(α1,α1)(α2,α1)...(αk,α1)(α1,α2)(α2,α2)...(αk,α2)............(α1,αk)(α2,αk)...(αk,αk)⎠⎟⎟⎞
Δ ( α 1 , α 2 , ⋯ , α k ) = ( α 1 α 1 T α 1 α 2 T . . . α 1 α k T α 2 α 1 T α 2 α 2 T . . . α 2 α k T . . . . . . . . . . . . α k α 1 T α k , α 2 T . . . α k , α k T ) \Delta (α_1,α_2,⋯,α_k)= \left( \begin{matrix} α_1α_1^T & α_1α_2^T & ... &α_1α_k^T\\ α_2α_1^T & α_2α_2^T & ... & α_2α_k^T\\ ... &...& ... &...\\ α_kα_1^T & α_k,α_2^T & ... & α_k,α_k^T \end{matrix} \right) Δ(α1,α2,⋯,αk)=⎝⎜⎜⎛α1α1Tα2α1T...αkα1Tα1α2Tα2α2T...αk,α2T............α1αkTα2αkT...αk,αkT⎠⎟⎟⎞
5.数据及代码下载地址
- GitHub的数据及代码下载地址为:GitHub的数据及代码下载链接(如果从GitHub下载代码,麻烦给小
Demo
一个Star
,您的支持是我最大的动力)