在线性SVM中，求解内积矩阵（也被称作Gram矩阵）
\[Q=\begin{bmatrix} \pmb x_i^T\pmb x_j \end{bmatrix}\]是一个必经的步骤，然后才可以进行后续的优化等步骤。当笔者在手动实现SVM时，就在求解上述$Q$矩阵时遇到困难。众所周知，Python在功能性强的同时舍弃了运行速度；在求得$Q$矩阵时，最基础的想法必然是通过二重循环实现：
# X是numpy矩阵，作为输入数据
Q = np.array([X[i] @ X[j] for j in range(len(X)) for i in range(len(X))])
我们来计算一下这样操作的时间，以sklearn自带数据集的digits数据集为例：
from sklearn.datasets import load_digits
import numpy as np
from time import time
X = load_digits().data
t = time()
for i in range(10):
Q = np.array([[X[i] @ X[j] for j in range(len(X))] for i in range(len(X))])
print((time() - t) / 10)
平均下来计算一次的时间是5.42秒，digits数据集是一个1797*64的数据集，这样小的数据集在循环计算下都需要5s附近，那么在大数据集下，这样的计算是及其耗时的。假如我们的数据是100000条digits数据集，也就是$10^{10}$次循环，我们就需要
的时间，我们无法接受。
重新回归计算过程，我们进行的其实是三重循环：
for i in range(1797):
for j in range(1797):
Q[i][j] = 0
for k in range(64):
Q[i][j] += X[i][k] * X[j][k]
这样的操作是不是很熟悉？没错，这就是矩阵乘法的伪代码，是一个$1797\times 64$矩阵和$64\times 1797$矩阵的乘积，显然就是$XX^T$，因此我们可以直接使用numpy中的矩阵乘积去实现，并计算运行时间：
t = time()
for i in range(10):
Q = X @ X.T
print((time() - t) / 10)
平均下来是0.0167秒，速度提升了300倍以上，而且形式简洁。
核方法通过核函数将输入空间上升到高维，从而解决一些线性模型无法解决的问题。理论上对于核函数可以分解：
\[K(\pmb x,\pmb y)=\phi(\pmb x)\phi(\pmb y)\]但实际上$\phi$函数往往不是显式的，即使显式的，作为一个由低维向高维投影的映射，其生成的向量也是长度极大的，会拖累计算速度。常用的核函数：
前面我们求解的就是线性核；从而，对于多项式核和Sigmoid核，只需要在线性核的基础上使用numpy的广播算法即可：
import numpy as np
gamma, theta, degree = 1, 0, 3
X = load_digits().data
linear = X @ X.T # 线性核
poly = np.power(gamma * linear + theta, degree) # 多项式核
sigmoid = np.tanh(gamma * linear + theta) # Sigmoid核
至于RBF核，计算关键是给定数据矩阵$X$求出$\Vert X_i-X_j\Vert_2^2$，按照之前的思路，我们希望能通过numpy数组的广播运算，求出这样的一个三维张量（不再是矩阵）：
也就是一个$N\times N\times d$的张量。直接相减会报错：
Y = X - X.T
# ValueError: operands could not be broadcast together with shapes (1797,64) (64,1797)
因此，我们将数据“升维”，然后利用广播算法相减：
Y = np.expand_dims(X, axis=1) - X
print(Y.shape) # 输出(1797, 1797, 64)
验证一下：
Y[1, 2] == X[1] - X[2] # 输出都为True，说明计算正确
然后利用numpy.linalg.norm计算每个$\Vert\pmb x_i-\pmb x_j\Vert_2^2$，形成一个$N\times N$矩阵，最后用广播算法算出RBF核函数：
Q = np.exp(-gamma * np.linalg.norm(np.expand_dims(X, axis=1) - X, axis=-1))
同样，我们可以像上面那样计算时间，时间比较：
| 核种类 | 时间(秒) |
|---|---|
| 线性核 | 0.01642 |
| 多项式核 | 0.13821 |
| Sigmoid核 | 0.04741 |
| RBF核 | 2.01980 |
可见RBF核计算代价是很大的。
再次考虑高斯核函数：
\[\begin{aligned} \exp(-\gamma\Vert\pmb x-\pmb y\Vert_2^2)&=\exp(-\gamma(\pmb x-\pmb y)^T(\pmb x-\pmb y))\\ &=\exp\big(-\gamma(\pmb x^T\pmb x+\pmb y^T\pmb y-2\pmb x^T\pmb y)\big) \end{aligned}\]因此，我们可以通过三次求内积矩阵来算出$\Vert\pmb x-\pmb y\Vert_2^2$：
gaussian_kenel = lambda x, y: np.exp(-gamma * ((x**2).sum(1).reshape(-1, 1) +
(y**2).sum(1) - 2 * x @ y.T))
比较两种方法的计算速度：
import numpy as np
from time import time
from sklearn.datasets import load_digits
gamma = 1
X = load_digits().data
kernel1 = lambda x, y: np.exp(-gamma * np.linalg.norm(
np.expand_dims(x, axis=1) - y, axis=-1)**2)
kernel2 = lambda x, y: np.exp(-gamma * ((x**2).sum(1).reshape(-1, 1) +
(y**2).sum(1) - 2 * x @ y.T))
t = time()
kernel1(X, X)
print(time() - t)
t = time()
kernel2(X, X)
print(time() - t)
| 计算方法 | 构造digits数据集核矩阵所需时间(秒) |
|---|---|
| 基于三维张量求解 | 2.27 |
| 基于代数求解 | 0.24 |
新方法大幅度加快了核函数矩阵的求解速度。在实现SVM时，我们将上面的算法投入了实践，并节省了时间。
