机器学习基础

本笔记是对 CMU Pratical Data Science Course 的机器学习相关部分的整理。

1. 一般化的机器学习

一般地，一个机器学习算法包含如下的核心要素：

• 输入(Inputs / Features) $x(i) \in \mathbb{R}^n, i=1, 2, \dots,m$。这是我们传入算法的信息。$x(i)$ 代表第 $i$ 个数据样本。它是一个 $n$ 维的实数向量。
• 输出 (Outputs) $y(i) \in Y$：这是我们希望算法预测的目标。$Y$ 代表所有可能输出的集合（输出空间）。输出可以是多种多样的，例如：
  • 回归 (Regression)：输出是连续的实数值（如用电量、房价）。
  • 分类 (Classification)：输出是离散的类别（如“是/否”、“猫/狗/鸟”）。
• 参数 (Parameters) $\theta \in \mathbb{R}^d$：$θ$ 是一个 $d$ 维向量，包含了模型中所有需要通过学习来确定的数值。我们的最终目的就是找到一组最优的 $θ$，使得模型的预测尽可能准确。
• 假设函数 (Hypothesis Function) $h_θ: R^n \rightarrow \hat{Y}$：这是一个从输入空间 $R^n$ 到预测空间 $\hat{Y}$ 的映射，就是我们的预测模型。给定一个输入 $x$，它会输出一个预测值 $hθ(x)$。
  • 模型的输出也由参数 $θ$ 决定，所以严格来说应该写作 $h(θ, x)$，但通常写作 $h_θ(x)$ 来强调它是一个从 $x$ 到 $y$ 的映射。
  • 预测空间 $\hat{Y}$ 不一定和真实输出空间 $Y$ 完全一样，但可以很容易地转换过去。例如在分类问题中，$h_θ(x)$ 可能输出一个0到1之间的概率，而真实的 $y$ 是0或1，这时两者就不同了。
• 损失函数 (Loss Function) $\ell: \hat{Y} × Y \rightarrow R_+$：这是一个衡量单次预测好坏的函数。它接收一个预测值 $h_θ(x)$ 和一个真实值 $y$，然后输出一个非负实数。

有了这些核心概念后，几乎所有的监督式机器学习算法都可以被归纳为以下的标准形式：

因此，要定义任何一个具体的机器学习算法，我们只需要回答以下三个问题：

1. 假设函数 $h_θ$ 的形式是什么？ (是线性模型、神经网络还是决策树？)
2. 损失函数 $\ell$ 是什么？ (是用平方误差、绝对值误差还是其他的？)
3. 我们如何解决这个最小化优化问题？ (是用梯度下降、解析解还是其他优化方法？)

下面我们从这些概念入手，讲解一些具体的方法。

2. 前置方法介绍

$a.$ 梯度下降

求解最小化优化问题的一个常用方法就是梯度下降。

$i.$ 简单例子

以下面这个假设函数为例：

我们计算待优化函数 $E(\theta)$ 相对 $\theta_1$ 和 $\theta_2$ 的梯度：

然后，我们就可以沿着梯度下降方向、以 $\alpha$ 步长更新这些参数、让 $E(\theta)$ 变小：

这里的 $E(\theta)=\sum_{i=1}^m(h_{\theta}(x^{(i)}-y^{(i)}))^2$ 叫作最小二乘估计。

更新 $\theta_1$，$\theta_2$ 简单代码如下：

theta = np.array([0., 0.])
alpha = 1.0
for t in range(100):
    theta[0] -= alpha/len(x) * 2 * sum((theta[0] * x_norm + theta[1] - y_norm)*x_norm)
    theta[1] -= alpha/len(x) * 2 * sum((theta[0] * x_norm + theta[1] - y_norm))

$ii.$ 最小二乘梯度下降方法整理

我们考虑 $E(\theta)$ 对 $\theta_j$ 的梯度：

这样，针对最小二乘法的梯度下降算法可以整理如下：

• 初始化：$\theta := 0$
• 迭代（直至收敛）：$For\;j\;=1,\dots,n$：

$b.$ 归一化

$i.$ 具体步骤

我们一般使用下面这个简单的操作，将 $\tilde{x}$ 映射到 $[0, 1]$ 区间中：

$\tilde(y)$ 同理。

$ii.$ 数据转换

我们在将数据归一化、进行梯度下降或者其他算法计算后，需要将它转换回原始数据。原始数据与归一化后数据的关系如下：

其中：

并且：

对于原始数据矩阵 $X$ 与结果向量 $y$，有：

其中新的系数为：

这一过程的简单代码实现如下：

def normalize_data(X, y, normalize_cols):
    """ Normalize y and specified columns of X in place. """
    # Normalize y in col
    min_X = X[:,normalize_cols].min(axis=0)
    max_X = X[:, normalize_cols].max(axis=0)
    min_y = y.min()
    max_y = y.max()

    X[:, normalize_cols] = (X[:,normalize_cols] - min_X) / (max_X - min_X)
    y[:] = (y - min_y) / (max_y - min_y)
    return (min_X, max_X, min_y, max_y)

def unnormalize_theta(theta, normalize_cols, ranges):
    # range tuple: (min_X, max_X, min_y, max_y)
    theta[normalize_cols] /= (ranges[1] - ranges[0])
    theta *= ranges[3] - ranges[2]
    theta[-1] += ranges[2] - theta[normalize_cols] @ ranges[0]

$c.$ 解析解

对于最小二乘法这个特殊形式，处理梯度下降方法，我们还可以采用解析解的方法求解最优状态。

解析解方法基于下面的基本原理：$E(\theta)$的最小值在梯度为0的地方取得。

先回到我们在最小二乘法的梯度下降公式：

这里和 $j$ 有关的项只有 $x_j^{(i)}$，于是我们整理一下这个式子的形式：

而根据线性代数中点积的定义：$\sum _{k=1} ^n \theta_kx_k=\theta^Tx=x^T\theta$，于是：

我们令 $\nabla E(\theta)=0$：

这个式子可以通过下面的代码计算：

def analytic_ls(X, y):
    m, n = X.shape
    A = np.zeros((n, n))
    b = np.zeros(n)
    for i in range(m):
        A += X[i].T @ X[i]
        b += X[i].dot(y)
    return np.linalg.solve(A, b)

$d.$ 矩阵记号下的梯度下降与解析解

我们将输入数据和期望结果用如下记号表示：

在这样的记号下，先前的最小二乘梯度计算：

可以重新写作：

这个公式十分简洁，同时解析解的形式也非常简洁：

两者的计算都只需要一行代码即可：

errors = X.T @ theta - y
grad = (2/m) * (X.T @ errors)

theta = np.linalg.solve(X.T @ X, X.T * y)

梯度计算式 $\nabla E(\theta) = \frac{2}{m}X^T(X\theta-y)$ 可以拆分成三个模块来理解：grad = (尺度调节器) * (责任分配器) @ (误差方向盘)。 $X\theta-y$ 计算了我们当前模型 $θ$ 对每一个数据点的预测值 $Xθ$ 与真实值 $y$ 之间的差距，这个值告诉我们误差的幅度和方向。

$X^T$ 则实现了将 $errors$ 误差进行分配的功能。由最初的梯度计算式，$θ_j$ 的总梯度，等于每一个样本的误差 $errors_i$，乘以该样本上 $θ_j$ 所对应的特征值 $X_{i,j}$，然后把这些乘积全部加起来。这就像一个加权投票的过程：$X$ 中的每个样本 $i$ 都在对 $θ_j$ 应该增加还是减少？”进行投票。投票的“票面价值”就是 $X_{i, j}$，投票的“权重”就是 $errors_i$

$\frac{2}{m}$ 是一个简单的“尺度调节器”，其中 $m$ 负责取平均值、让梯度更加稳定；而 2 只是一个求导过程产生的参数，可以被学习率 $\alpha$ 吸收。

$e.$ 线性回归简单实现

根据上面的内容，我们仿照 scikit-learn 的设计，实现一个简单的线性回归类：

class MyLinearRegression:
    def __init__(self, method='gd', learning_rate=0.1, n_iters=1000):
        self.method = method
        self.learning_rate = learning_rate
        self.n_iters = n_iters
        self._theta = None
    
    def fit(self, X, y):
        # We add bias item to the matrix 
        m, n = X.shape
        X_b = np.c_[np.ones((m, 1)), X]

        if self.method == 'gd':
            theta = np.zeros(0)
            for i in range(self.n_iters):
                errors = X_b @ theta - y
                grad = (2/m) * (X_b.T @ errors)
                theta = theta - self.learning_rate * grad
            self._theta = theta

        elif self.method == 'normal':
            A = X_b.T @ X_b
            b = X_b.T @ y
            self.theta_ = np.linalg.solve(A, b)
        else:
            raise ValueError("method must be gd or normal")
        
        return self 
    
    def predict(self, X):
        if self.theta_ is None:
            raise RuntimeError("Please train the model before prediction")
        
        X_b = np.c_[np.ones((X.shape[0], 1)), X]
        return X_b @ self.theta_
    
    @property
    def intercept_(self):
        return self.theta_[0] if self.theta_ is not None else None
    
    @property
    def coef_(self):
        return self.theta_[1:] if self.theta_ is not None else None

3. 线性分类问题

$a.$ 相关概念

与前面讨论的回归问题相比，分类问题有如下特点：

• 输出空间 $Y$ 是有限的、离散的集合，而不是连续的区间。
• 预测空间 $\hat{Y}$ 仍然是连续的实数 $\mathbb{R}$。在中间过程中，我们的假设模型 $h_{\theta}$ 会输出一个实数（如5.2，-1.7）,其中：
  • 值的符号 (sign) 决定了最终的预测类别：如果 hθ(x) 是正数，就预测为 +1；如果是负数，就预测为 -1。
  • 值的绝对值大小 (magnitude) 代表了预测的把握有多大。

通过让我们的假设模型返回一个连续的实数结果，我们可以得到更为平滑的结果、更好地进行梯度下降等算法。因为 $\lbrace +1, -1\rbrace$ 这样的输出是非连续、不平滑的。

$b.$ 损失函数设计

分类问题常用的损失函数如下：

• 逻辑损失 (Logistic Loss) $\log(1 + \exp(-h_θ(x) \cdot y))$

  • 当预测非常正确时，$h_θ(x)\cdot y$ 是很大的正数），$\exp$ 趋近于0，总损失趋近于 $\log(1)=0$。
  • 当预测非常错误时，$h_θ(x)\cdot y$ 是很大的负数），损失近似线性增长。

• 合页损失 (Hinge Loss) $\max \lbrace 1 - h_θ(x) \cdot y, 0 \rbrace$

  • 只要 $h_θ(x) \cdot y \geq 1$（即不仅分类正确，而且置信度超过一个“边界”1），损失就精确地为 0。
  • 如果 $h_θ(x) \cdot y \lt 1$，损失会线性增加。

• 指数损失 (Exponential Loss) $\exp(-h_θ(x) \cdot y)$

  • 当预测正确时，损失趋近于 0。
  • 当预测错误时，损失会指数级地爆炸式增长。

有了这些损失函数，我们的分类问题处理框架就和回归问题的类似了：

$c.$ 支持向量机

支持向量机 (Support vector machines, SVMs) 算法使用合页损失。SVM的假设函数有线性和核性两种。线性 SVM 使用如下的假设函数：

核性 SVM 使用如下的假设函数：

添加了正则化参数的完整表达式如下：

我们推导 SVM 梯度的计算。假设函数部分：

正则化部分：

于是：

矩阵形式的梯度表达式如下：

其中 $Y=diag(y)$。

这里的 $Y X\theta \le 1$ 可以理解为“只对支持向量进行计算，其他样本的梯度贡献为0”，这也是“支持向量机”算法名称的来源。

根据这个矩阵形式梯度计算公式，我们可以写出如下代码：

def svm_dg(X, y, alpha, iters, lam):
    m, n = X.shape
    theta = np.zeros(n)
    YX = X * y[:, None]
    loss, err = np.zeros(iters), np.zeros(iters)

    for i in range(iters):
        hy = YX @ theta
        loss[i] = np.maximum(1-hy, 0).mean()
        err[i] = (hy <= 0).mean()
        theta = theta - alpha * (-YX.T @ (YX * theta <= 1) / m + 2 * lam * theta)

    return theta, loss, err

这里的 YX.T 可以通过维度匹配判断需要转置。

$d.$ 逻辑回归

逻辑回归 (Logistic Regression) 算法使用逻辑损失函数：

我们推导逻辑回归梯度的计算：

因此：

矩阵形式的梯度表达如下：

根据这个矩阵形式梯度计算公式，我们可以写出如下代码：

def logistic_reg_gd(X, y, alpha, iters):
    m, n = X.shape
    theta = np.zeros(n)
    YX = X * y[:, None]
    loss, err = np.zeros(iters), np.zeros(iters)

    for i in range(iters):
        hy = YX @ theta
        loss[t] = np.log(1+np.exp(-hy)).mean()
        err[t] = (hy <= 0).mean()
        theta -= alpha * -YX.T @ (1/(1+np.exp(YX @ theta)))/m
    return theta, loss, err

4. 多种类分类问题

前面我们讨论的都是二元分类问题，但是现实生活中有很多问题是需要多重分类的。多重分类的输出空间是一个包含 $k$ 个离散值的集合：$Y=\lbrace 1, 2, \dots,k \rbrace$。处理多分类问题一般有下面两种方法：

• “一对多” (One-versus-All, OvA)：把一个 $k$ 分类问题，拆解成 $k$ 个独立的二元分类问题。具体做法如下：
  1. 为每一个类别 $i$，我们都训练一个专属的二元分类器 $h_θi$。这个分类器 $h_θi$ 的任务是回答一个“是/否”问题：“这个样本是第 $i$ 类，还是不是第 $i$ 类？”
  2. 在为分类器 $h_θi$ 准备训练数据时，我们将所有原始标签为 $i$ 的样本，其新标签设为 +1，其他样本的新标签设为 -1。
  3. 重复这个过程 $k$ 次，我们就得到了 $k$ 个训练好的二元分类器。

当来一个新的样本 $x$ 时，我们让它分别通过所有 $k$ 个分类器，得到 $k$ 个代表“置信度”的实数值输出：$h_{θ1}(x), h_{θ2}(x), ..., h_{θk}(x)$。由于可能有多个分类器返回整数，因此我们比较它们的绝对大小，选择给出了最高置信度分数的分类器所对应的类别：

• 原生多分类 (Native Multiclass)：我们可以不把 $k$ 个分类器看作是独立的，而是把它们看作一个统一的、更复杂的分类器。这个统一的分类器的假设函数 $h_θ$，其输出不再是一个单独的实数，而是一个 $k$ 维的向量 $R^k$。$k$ 维向量的第 $i$ 个元素 $h_θ(x)^i$，就代表了模型认为样本属于第 $i$ 类的“相对置信度”。

有了 $k$ 维的输出，我们就可以定义一个“原生”的损失函数 $\ell$。这个函数接收一个 $k$ 维的预测向量 $h_θ(x)$ 和一个真实的类别标签 $y$，然后计算出一个损失值。常见的损失函数有 softmax 函数

5. 泛化与过拟合

$a.$ 相关概念

在前面内容中，我们提到在训练机器学习模型的一个任务：

然而，这只是手段，不是最终目的。因为训练集里的数据是已知的，我们已经知道正确答案 $y(i)$ 了。我们真正关心的是泛化能力 (Generalization)：当一个全新的、模型从未见过的数据出现时，我们的模型能否依然做出准确的预测。

由此，我们引入泛化误差 (Generalization Error) 。它指的是模型在新数据上的预测误差。我们真正需要的是最小化泛化误差，而不是训练误差。

而一个模型如果在训练集上表现极好（训练误差很低），但在新数据上表现很差（泛化误差很高），它就是过拟合 (Overfitting) 的。

$b.$ 交叉验证

交叉验证是找到那个泛化误差最低的最佳点的常用手段。

最简单的交叉验证方法如下：

1. 将原始数据集分成两部分，例如：70% 作为训练集 (Training Set)，30% 作为验证集 (Validation Set / Hold-out Set)。
2. 我们只用训练集来训练模型（即调整参数 $θ$）。
3. 训练完成后，用验证集来测试模型，计算其误差（称为验证误差）。

验证集和测试集在概念上是不同的。验证集用于选择模型和调整超参数，而测试集是在所有选择和调整都完成后，对最终模型的一次性最终评估，评估后就不能再返回去修改模型了。

简单的70/30划分有一个缺点：结果可能受到单次随机划分的影响。如果某次划分碰巧把一些很“奇怪”的数据点（离群点）都分到了验证集，那么评估结果可能就不太准。为了解决这个问题，我们可以使用K-折交叉验证，它的步骤如下：

1. 将整个数据集平均分成 $k$ 份（称为 k fold）。
2. 进行 $k$ 次训练和验证，每次：
   • 选其中 1 份作为验证集。
   • 剩下的 $k-1$ 份作为训练集。
3. 将 $k$ 次验证得到的误差取平均值，作为最终的评估结果。

6. 正则化

$a.$ 相关概念

机器学习模型有一个重要的参数：模型参数（系数）的大小 (magnitude of the model parameters)。一个模型如果依赖巨大的、精巧抵消的系数，说明它非常不稳定和复杂。它对输入数据的微小变化会非常敏感。一个更简单、更鲁棒的模型，其系数通常会比较小，曲线也更平滑。

正则化正是一种在机器学习模型的优化过程中，通过在损失函数中加入一个惩罚项 (penalty term) 来限制模型参数大小的技术。

$b.$ 最小二乘中的正则化

我们将 $\ell_2$ 正则化应用到最小二乘模型中:

在前面的内容中，我们已经知道平均损失部分的梯度：

而正则化部分的梯度：

于是：

我们可以把 $m$ 并入正则化参数 $\lambda$ 中：

这个结果体现了正则化的另一个好处：在标准的最小二乘法中，如果 $X^TX$ 不可逆（比如特征之间高度相关），就无法求解。但加入正则化项后，由于我们加上了一个 $λI$ 项（即在 $X^TX$ 的对角线上都加上一个小的正数 $\lambda$），这使得矩阵几乎总是可逆的，从而求解过程更加稳定。

这个正则化项可以简单地迁移到梯度下降中，我们只需要将参数更新规则改为下面的形式：

6. 非线性特征

$a.$ 特征映射

特征映射是将非线性问题转化成线性问题的一个常用方法。特征映射将 $n$ 维原始输入映射为例 $k$ 维的特征向量：

这样，我们就可以转而考虑一个线性的多项式函数（对 $k$ 阶向量而言）作为假设函数，它的参数 $\theta\in\mathbb{R}^k$:

特征映射的神奇之处在于：我们不需要发明新的、复杂的非线性算法，我们只需要在应用简单线性算法之前，对数据进行一次“升维和变换”。

$b.$ 径向基函数

$i.$ 定义与特征

径向基函数的形式如下：

其中 $\mu^{(1)},\dots,\mu^{k-1}$ 称作均值，$\sigma$称作带宽。

与多项式这样的全局特征 (global features) 不同，RBF 是局部特征 (local features)。RBF 的每个特征只在输入空间的一小块区域内有显著的非零值，在其他地方都接近于零。

这些参数对特征的影响如下：

• $μ(j)$：是第 $j$ 个基函数的中心 (center)。它决定了这个“钟形”的峰值在哪里。
• $σ$: 它决定了这个“钟形”的胖瘦。$σ$ 越大，钟形越胖，特征影响的范围越广；$σ$ 越小，钟形越瘦，特征影响的范围越窄。
• $\exp(-...)$: 确保当 $x$ 离中心 $μ(j)$ 越远，$φ_j(x)$ 的值就越快地衰减到 0。当 $x = μ(j)$ 时，$φ_j(x)$ 取最大值 1。

和多项式一样，RBF 也会过拟合。如果我们使用的 RBF 特征数量过多（即中心 $μ$ 的数量太多），模型就会变得过于灵活，试图去拟合训练数据中的每一个点，导致函数曲线变得非常“扭曲”和“振荡”。

最终 RBF 会得到一个如下的特征矩阵：

然后我们就可以使用这个矩阵作为特征输入，进行梯度计算等操作了。

使用 RBF 作为特征工程的简单实现如下：

def rbf_feat(x, mu, sig):
    return np.hstack([np.exp(-(x[:, None] - mu)**2 / (2 * sig ** 2)), np.ones(len(x), 1)])

def train_rbf(x, y, n_rbf):
    min_x, max_x = x.min(), x.max()
    sig = (max_x - min_x)/(n_rbf-1)
    mu = np.linspace(min_x, max_x, n_rbf-1)
    
    Phi = rbf_feat(x, mu, sig)
    theta = np.linalg.solve(Phi.T @ Phi, Phi.T @ y)
    return theta, mu, sig

$ii.$ 超参数

RBF 的主要超参数如下：

• 中心点的数量和位置 $μ$：这是最直接控制模型复杂度的参数。中心点越多，模型越复杂。
  • 在一维数据中，可以简单地将它们均匀分布。但在高维数据中，如何选择中心点本身就是一个问题（例如，可以随机选一些数据点作为中心，或使用聚类算法来确定中心）。
• 正则化参数 $λ$
• 带宽参数 $σ$：带宽也控制着模型的复杂度。当 $\sigma$ 较大时，钟形曲线很胖，各个基函数之间重叠很多。最终的函数会更平滑，模型更简单。增大带宽 $σ$ 和增大正则化参数 $λ$ 有类似的效果，都是在简化模型，使其更平滑。

7. 多维度中的非线性特征

$a.$ 多维多项式

当原始输入 $x$ 是一个 $n$ 维向量时，如果需要构建它的多项式特征，除了每个维度自身的高次项 $x_1^2$, $x_2^2$, $x_3^3$ 外，我们还必须包含不同维度之间的交叉项 $x_1x_2$, $x_1^2x_2$, $x_1x_2^2$ 等。

Python 的 itertools 模块的 combinations_with_replacement 方法可以从一个列表中有放回地选出指定数量的元素。利用这个方法可以轻松生成 $x_1,\dots,x_n$ 的各种组合：

• 如果我们需要生成所有阶为 $n$ 的项：

x = np.array(["x1","x2","x3"])
list(combinations_with_replacement(x,3))

# output:
# [('x1', 'x1', 'x1'),
#  ('x1', 'x1', 'x2'),
#  ('x1', 'x1', 'x3'),
#  ('x1', 'x2', 'x2'),
#  ('x1', 'x2', 'x3'),
#  ('x1', 'x3', 'x3'),
#  ('x2', 'x2', 'x2'),
#  ('x2', 'x2', 'x3'),
#  ('x2', 'x3', 'x3'),
#  ('x3', 'x3', 'x3')]

• 如果我们需要生成最高阶为 $n$ 的项，我们可以在数组中添加一个常数项 1：

x = np.array(["x1","x2","x3","1"]) # just an example
list(combinations_with_replacement(x,3))

# output:
# [('x1', 'x1', 'x1'),
#  ('x1', 'x1', 'x2'),
#  ('x1', 'x1', 'x3'),
#  ('x1', 'x1', '1'),
#  ('x1', 'x2', 'x2'),
#  ('x1', 'x2', 'x3'),
#  ('x1', 'x2', '1'),
#  ('x1', 'x3', 'x3'),
#  ('x1', 'x3', '1'),
#  ('x1', '1', '1'),
#  ('x2', 'x2', 'x2'),
#  ('x2', 'x2', 'x3'),
#  ('x2', 'x2', '1'),
#  ('x2', 'x3', 'x3'),
#  ('x2', 'x3', '1'),
#  ('x2', '1', '1'),
#  ('x3', 'x3', 'x3'),
#  ('x3', 'x3', '1'),
#  ('x3', '1', '1'),
#  ('1', '1', '1')]

在输入矩阵 $X$ 中，我们可以通过下面这个简单函数计算矩阵元素所有多项式乘积：

def poly_feat(X,d):
    m = X.shape[0]
    X_ = np.hstack([X, np.ones(m, 1)]).T
    return np.array([np.array(a).prod(axis=0) for a in combinations_with_replacement(X_,d)]).T

在数学计算中，一个维度为 $n$ 的 $d$ 阶多项式具有的特征数量为 $C(n+d), d$。

$b.$ 多维 RBF

为了将 RBF 推广到多维，只需要将一维的距离 $(x - μ)^2$ 替换为多维空间中的平方欧氏距离 $||x - μ||^2$：

直接计算平方欧式距离的效率较低，我们可以通过矩阵向量运算进行优化，对单个向量我们有：

由于我们最终需要得到维度为 $m \times k$ 的矩阵，因此 $x^{(i)T} \mu^{(j)}$ 转换成矩阵形式为 $XM^T$。利用 numpy 的广播机制，我们可以写出如下的函数：

def sqdist(X, M):
    return (X**2).sum(axis=1)[:,None] + (M**2).sum(axis=1) - 2*X @ M.T

同样地，多维 RBF 特征计算如下：

def rbf_feat(X, mu, sig):
    return np.exp(-sqdist(X,mu) / (2*sig**2))

在数学计算中，$n$ 维 $d$ 中心的 RBFs 具有的特征数量为 $d^n$。

$c.$ 多项式与 RBF 多维特征的实用准则

多维多项式特征的一些基本准则如下：

• $C(n+d, d)$ 这个公式告诉我们，只有当维度 $n$ 和阶数 $d$ 同时很大时，特征数量才会爆炸。因此：
  • 如果原始维度 $n$ 很高，那么只使用很低的阶数 $d$（通常 $d=2$，即只生成二次项，这已经很常见且有效了）。
  • 如果确实需要使用很高的阶数 $d$，那么必须保证原始维度 $n$ 非常低。
• 处理二元特征 (Binary Features)：如果一个特征 $x_j$ 是二元的（即取值只有 0 或 1），那么 $x_j^2 = x_j$, $x_j^3 = x_j$，这意味着对这个特征自身取更高次幂是毫无意义的，因为它们是完全冗余的。但是，它与其他特征的交叉项$x_j \cdot x_{other}$ 仍然是有意义的。
• 选择性地创建交叉项：不需要创建所有可能的交叉项，而是可以根据对问题的理解，只创建那些有意义的。

RBFs 的一些基本准则如下：

• $\mu$ 应该从数据本身中选取。真实世界的高维数据几乎从不均匀地分布在整个空间里，而是倾向于聚集在某个低维的流形 (manifold)上。我们应该让中心点 $\mu$ 分布在这些有实际数据的地方。
  • 一个简单的方法是从已有 $m$ 个数据点中，随机抽取 $k$ 个点作为 RBF 的中心 $\mu$。更优的方法是使用 k-means 聚类算法：对数据集运行 k-means，找到 k 个簇的质心 (centroids)，然后用这 k 个质心作为 RBF 的中心点。这是一个非常常用且有效的策略，因为它找到的点既能代表数据的高密度区域，又相对分散，能很好地覆盖住数据所在的流形。
• $\sigma$ 可以选取为数据的中位数：$σ = \text{median}(||x^{(i)} - μ^{(j)}|| i, j = 1, 2, \dots,k-1)$。中位数是数据的“典型值”的一个非常稳健的估计量，它不受少数极大或极小距离值的影响。

8. 核方法

核方法的巧妙之处在于：它允许我们在一个极高维度（甚至是无限维度）的特征空间中进行计算，而完全不需要显式地创建或存储这些高维特征向量。

$a.$ $\theta$ 的新表达形式

首先，我们提出表示定理：对于一个线性模型 $h(x) = θ^Tφ(x)$，其最优解 $\theta^\star$ 总可以表示为所有训练数据点特征向量的线性组合：

对这个定理的直观理解如下：最优的权重向量 $θ$ 一定存在于由训练数据的特征向量 $φ(x^{(i)})$ 所张成的空间中。任何超出这个空间的分量都是无用的，因为它与所有训练数据都正交，对损失函数没有贡献。

这个结论意味着，我们不再需要直接求解 $k$ 维的 $θ$，而是可以转而求解 $m$ 个权重 $α$（$m$ 是训练样本的数量）。

既然 $θ$ 可以用 $α$ 来表示，我们也可以把预测函数 $h(x)$ 也重写成只与 $α$ 相关的形式：

这个式子展示了一个重要的事实：对一个新点 $x`$ 的预测只依赖于它的特征向量 $φ(x)$ 与所有训练点的特征向量 $φ(x(i))$ 之间的点积。我们根本不需要知道 $φ(x)$ 的具体长什么样，我们只需要知道如何计算它与其他向量的内积。而这里就是核计算发挥作用的地方。

$b.$ 点积的高效计算

在核方法中，我们尝试找到一个函数 $K(x, z)$，它能直接在原始的、低维的输入空间中计算出高维特征向量的内积 $φ(x)^Tφ(z)$，从而完全绕过高维特征向量的显式计算。函数 $K(x,z)$ 就是核函数。

$i.$ 多项式核

对于一个 $d$ 阶多项式，我们有：

对应的核函数 $K: \mathbb{R}^n \times \mathbb{R}^n \rightarrow \mathbb{R}$：

$ii.$ RBF 核

RBF 特征的核函数如下：

RBF 核的简单实现如下：

def gaussian_kernel(X1, X2, sig):
    sqd = sqdist(X1, X2)
    return np.exp(-sqd/(2 * sig ** 2))

我们进一步讨论这个式子展现出的 RBF 的无限维度的表达能力。我们考虑下面的指数核：

这个指数核对应的特征向量 $\varphi(x)$、$\varphi(z)$ 可以通过泰勒展开得到：

这正是向量点积的形式，于是我们的特征向量可以写为：

这个特征向量 $φ(x)$ 包含 x 的所有次幂（从0次、1次、2次...一直到无穷），所以它是一个无限维的向量。

通过这个无限维的向量，我们获得了无限维度的强大能力，却付出了极低的计算代价。我们完全避免了“维度灾难”，因为我们从始至终都没有真正创建那个无限维的向量。这正是核方法最精妙、最强大的地方。