（1）Python之向量（Vector）距离矩阵计算

1. 引言

        距离矩阵是一个包含一组点两两距离的矩阵（即 二维数组）。因此给定$N$个欧几里得空间中的点，其距离矩阵就是一个非负实数作为元素的$N\times N$的对称矩阵。在机器学习中距离矩阵都计算非常常见（只要涉及距离计算，基本都需要计算距离矩阵），在本篇博客中就来记录一下如何使用Python都科学计算包numpy计算向量都距离矩阵。本篇博客讲解以行向量的欧氏距离为例讲解，但是同时给出了列向量的代码，距离矩阵的数学表达为：

给定$n\times m$阶矩阵$X$，满足$X=\left[\begin{array}{c}{x}_1\\ x_2\\ ...\\ x_n\end{array}\right]$，里第$i$行向量是$m$维向量**（注意这里是行向量）**。使得：

$$D_{ij}=\Vert x_i-x_j\Vert$$

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这个方法，这个方法从数学上讲，其计算公式是：

$$\Vert x_i-x_j\Vert =\sqrt{(x_i-x_j)(x_i-x_j{)}^T}$$
但是从另一方面来讲，我们可以先求点积运算，然后在进行求根运算：

$$D_{ij}=\sqrt{(x_i-x_j)(x_i-x_j{)}^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$(针对列向量，如果是行向量就是$m^2-m$)次，并且每次进行了m(针对列向量，如果是行向量就是$n$)次乘积运算和n-1(针对列向量，如果是行向量就是$m^2-m$)次加法运算。尽管numpy底层可能对点积运算做了优化，但这里还是存在可能进行进一步优化。请看下面的数学推导（行向量）：

$$D_{ij}=\sqrt{(x_i-x_j)(x_i-x_j{)}^T}=\sqrt{x_i{x}_i^T-2x_i{x}_j^T+x_j{x}_j^T}$$

这里$x_i{x}_i^T,x_j{x}_j^T,x_i{x}_j^T$属于格拉姆矩阵中的元素，可以通过在循环外计算矩阵，在循环内直接引用元素值即可，从而在循环内我们只需要做两次加（减）法运算：

$$D_{ij}=\sqrt{G_{ii}-2G_{ij}+G_{jj}}$$

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

# 行向量
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=\sqrt{H+K-2G}$，与公式$D_{ij}=\sqrt{G_{ii}-2G_{ij}+G_{jj}}$进行对比，有：

$$H_{ij}=G_{ii},K_{ij}=G_{jj}$$
这里H中第i行的每一个元素，取值都为$G_{ii}$，也就是H的每一列，都对应着格拉姆矩阵的对角阵，因此，我们可以用下面的代码来计算H(n是向量的个数，无论行向量还是列向量)（帮助理解：将G都对角线元素展开成一个行向量，那么$H_{ij}$的每一个元素只和i有关，即相应G位置向上找，所以将行向量数值排列）：
$$H=np.tile(np.diag(G),(n,1))\text{\hspace{1em}(将Y轴复制n倍，将X轴复制1倍)}$$
此外，由于$K=H^T$,所以最终距离矩阵可以计算为
$$D=\sqrt{H+H^T-2G}\text{\hspace{1em}(开根代表对矩阵每一个元素求根)}$$

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

# 行向量
#[[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\times 3072$阶矩阵。这里5000代表5000幅带标签的图，3072是其各像素在RGB三个通道下的取值数。给定测试集矩阵B为$500\times 3072$阶矩阵。求矩阵B的各行与矩阵A的各行的距离（即两幅图的差异）矩阵，这个矩阵是一个$5000\times 500$的矩阵。

更一般地，这个问题可以描述如下：

给定矩阵A为$m\times k$阶矩阵，矩阵B为$n\times k$阶矩阵，求矩阵B的任意行向量与矩阵A的任意行向量的距离矩阵$D_{m\times n}$。这个矩阵的数学表达式为($a,b$均为行向量）：
$$D_{ij}=\sqrt{(a_i-b_j)(a_i-b_j{)}^T}=\sqrt{a_i{a}_i^T-a_i{b}_j^T-b_j{a}_i^T+b_j{b}_j^T}$$

为方便讨论，我们将上述各项分别记为$H,M,N,K$，即：
$$D_{m\times n}=\sqrt{H_{m\times m}-M_{m\times n}-N_{n\times m}+K_{n\times n}}\text{\hspace{1em}(开根代表对矩阵每一个元素求根)}$$

显然上述公式是无法进行运算的，因为除了M与D外，其它矩阵的秩各不相同。所以我们要回到前一个数学表达式上。

1. $H$对$D$的贡献是对于$D$的每一行，都加上$a_i{a}_i^T$
2. $K$对$D$的贡献是对于$D$的每一列，都加上$b_j{b}_j^T$
3. $M$和$N$互为转置矩阵。即对$D_{ij}$,要减去矩阵$N_{j,i}$元素，而这个元素就是$M_{ij}$
4. 由前三个总结$D_{ij}$为$a_i{a}_i^T-a_i{b}_j^T-b_j{a}_i^T+b_j{b}_j^T$求根

代码如下（行向量）：

# 行向量: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,\cdots ,v_n$ 是内积空间的一组行(列)向量，Gram矩阵定义为： $G_{ij}=⟨v_i,v_j⟩=v_i{v}_j^T$，显然其是对称矩阵。
• 其实对于一个$X_{N\cdot d}$（行向量N 个样本，d 个属性）的样本矩阵而言，$X\cdot X^T$ 即为 Gram 矩阵。如果是列向量$X^T\cdot X$ 为 Gram 矩阵(很重要，这里需要用到的)
• 如果 $v_1,v_2,\cdots ,v_n$分别是随机向量，则 Gram 矩阵是协方差矩阵；
• 欧式空间中向量$v_1,v_2,\cdots ,v_n$的Gram矩阵一定是半正定矩阵，是正定矩阵的充要条件是$v_1,v_2,\cdots ,v_n$线性无关。

具体形式为：$n$维欧式空间中任意$k(k\le n)$个向量${\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_k$的内积所组成的矩阵
$$\Delta \left({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_k\right)=\left(\begin{array}{cccc}({\alpha }_1,{\alpha }_1)& ({\alpha }_1,{\alpha }_2)& ...& ({\alpha }_1,{\alpha }_k)\\ ({\alpha }_2,{\alpha }_1)& ({\alpha }_2,{\alpha }_2)& ...& ({\alpha }_2,{\alpha }_k)\\ ...& ...& ...& ...\\ ({\alpha }_k,{\alpha }_1)& ({\alpha }_k,{\alpha }_2)& ...& ({\alpha }_k,{\alpha }_k)\end{array}\right)$$

$$\Delta ({\alpha }_1,{\alpha }_2,\cdots ,{\alpha }_k)=\left(\begin{array}{cccc}{\alpha }_1{\alpha }_1^T& {\alpha }_1{\alpha }_2^T& ...& {\alpha }_1{\alpha }_k^T\\ {\alpha }_2{\alpha }_1^T& {\alpha }_2{\alpha }_2^T& ...& {\alpha }_2{\alpha }_k^T\\ ...& ...& ...& ...\\ {\alpha }_k{\alpha }_1^T& {\alpha }_k,{\alpha }_2^T& ...& {\alpha }_k,{\alpha }_k^T\end{array}\right)$$

5.数据及代码下载地址

• GitHub的数据及代码下载地址为：GitHub的数据及代码下载链接（如果从GitHub下载代码，麻烦给小Demo一个Star，您的支持是我最大的动力）