[Machine Learning] Feature Selection

---

Subset Selection

Forward Stepwise Selection

1. 初始化一个不包含任何预测变量的模型，这被称为null model。

2. 对于 $k=0,1,2,...,p-1$（这里 $p$ 是总预测变量的数量），执行以下步骤：
   

   • 考虑所有的 $p-k$ 个模型，这些模型将 $M_k$ 中的预测变量增加一个额外的预测变量。
   • 选择这些 $p-k$ 个模型中最好的一个，将其指定为 $M_{k+1}$。这里的“最好”是根据一些指标（例如，残差平方和（RSS）或分类错误）来衡量的。
     

3. 在所有的 $M_0,M_1,\dots ,M_p$ 模型中选择单个最好的模型。这可以通过交叉验证预测错误，残差平方和等进行。

Forward Stepwise Selection并不能保证在所有的 $2^p$ 个包含 $p$ 个预测变量子集的模型中找到最好的模型。这是因为在每一步中，它只选择在当前模型中增加一个预测变量能够最大程度地提升模型拟合度的模型，而不是考虑所有可能的预测变量子集。

Backward Stepwise Selection

1. 初始化一个包含所有预测变量的模型，这被称为full model。

2. 对于 $k=p,p-1,...,1$（这里 $p$ 是总预测变量的数量），执行以下步骤：

   • 考虑所有的 $k$ 个模型，这些模型包含 $M_k$ 中的所有预测变量，但删除了一个，总共有 $k-1$ 个预测变量。
   • 选择这些 $k$ 个模型中最好的一个，将其指定为 $M_{k-1}$。这里的“最好”是根据一些指标（例如，残差平方和（RSS）或分类错误）来衡量的。

3. 在所有的 $M_0,M_1,\dots ,M_p$ 模型中选择单个最好的模型。这可以通过交叉验证预测错误，残差平方和等进行。

与Forward Stepwise Selection类似，Backward Stepwise Selection只需要搜索 $1+p(p+1)/2$ 个模型，所以在 $p$ 太大以至于不能应用最佳子集选择的情况下，也可以应用Backward Stepwise Selection。

注意：对于一些模型（如线性回归），Backward Stepwise Selection要求样本数量 $n$ 大于特征数量 $p$ （以便能够拟合full model）。

尽管如此，Backward Stepwise Selection并不能保证在所有的 $2^p$ 个包含 $p$ 个预测变量子集的模型中找到最好的模型。这是因为在每一步中，它只选择在当前模型中删除一个预测变量能够最大程度地提升模型拟合度的模型，而不是考虑所有可能的预测变量子集。

Estimating Test Error

Mallow的$C_p$

Mallow的$C_p$是一种在多元线性回归模型中选择最优模型的方法。Mallow的$C_p$统计量以预测残差平方和和模型中参数数量为基础，这使得我们不仅可以评价模型的拟合优度，而且还可以对模型复杂度进行考虑。

Mallow的$C_p$的公式如下：

$$C_p=\frac{1}{n}(RSS+2d{\hat{\sigma }}^2)$$

其中：
$n$ 是样本数；
$RSS$ 是残差平方和，定义为 ${\sum }_{i=1}^n(y_i-{\hat{y}}_i{)}^2$；
$d$ 是模型中的参数数量；
${\hat{\sigma }}^2$ 是对方差的估计，通常由全模型（即包含所有预测变量的模型）的残差平方和除以对应的自由度得到。

在这个公式中，第一个部分 $\frac{1}{n}RSS$ 表示模型拟合的优度，第二个部分 $2d{\hat{\sigma }}^2$ 用来对模型复杂度进行惩罚。模型复杂度越高（即参数数量 $d$ 越多），模型可能就越倾向于过拟合，所以需要通过复杂度惩罚来对这种可能性进行调整。

使用Mallow的$C_p$选择模型时，通常会选择$C_p$最小的模型，因为这表示模型既有良好的拟合度，又没有引入过多的复杂性。

贝叶斯信息准则 (Bayesian Information Criterion, BIC)

$$BIC=n\ln (\frac{RSS}{n})+k\ln (n)$$

其中：

$n$ 是样本数量
$RSS$ 是模型的残差平方和
$k$ 是模型中的参数数量

在实践中，我们通常选择BIC最小的模型作为最优模型。这是因为BIC的最小化相当于最大化模型的后验概率，假设模型参数的先验分布是均匀的。

Adjusted $R^2$

$$\overline{R^2}=1-[(1-R^2)\frac{n-1}{n-p-1}]$$

其中：
$R^2$ 是未调整的决定系数；
$n$ 是样本数量；
$p$ 是模型中的参数数量（不包括截距项）。

在实践中，我们通常选择调整的$R^2$最大的模型作为最优模型。

Shrinkage

Ridge回归 (Tikhonov Regularization)

在标准的线性回归中，我们的目标是最小化RSS：

$$RSS=\sum _{i=1}^n(Y_i-{\beta }_0-\sum _{j=1}^p{\beta }_j{X}_{ij}{)}^2=||Y-X\beta |{|}_2^2$$

这里，$X$ 是我们的特征矩阵，$y$ 是目标向量，$\beta$ 是我们希望估计的系数向量。可以看到，如果我们将 $X$ 中的一个特征列 $X_j$ 乘以一个常数 $c$ ，那么，求解 $\beta$ 时，对应的 ${\beta }_j$ 就会除以同样的常数 $c$。因此，我们说最小二乘法的系数估计是尺度不变的。

对于岭回归，我们的目标函数是最小化残差平方和与系数向量的 $L_2$ 范数的和：

$$RSS+\lambda \sum _{j=1}^p{{\beta }_j}^2=||Y-X\beta |{|}_2^2+\lambda ||\beta |{|}_2^2$$

这里，$\lambda$ 是一个非负的调整参数，用于控制正则化的强度。在这种情况下，如果我们将 $X$ 中的一个特征列 $X_j$ 乘以一个常数 $c$ ，那么，求解 $\beta$ 时，对应的 ${\beta }_j$ 的变化就会受到 $c$ 的影响，因为 $\lambda ||\beta |{|}^2$ 部分也受到 ${\beta }_j$ 的影响，导致对应的系数估计会有实质性的变化。因此，我们说岭回归的系数估计是尺度变化的。

这也就解释了为什么在使用岭回归之前，我们通常需要对特征进行标准化（即，将特征的平均值变为0，标准偏差变为1），这样所有特征都在同一尺度上，不会因为尺度的差异而影响系数的估计。特征标准化如下：

$${\tilde{X}}_{ij}=\frac{1}{\sqrt{\frac{1}{n}{\sum }_{i=1}^n(X_{ij}-{\bar{X}}_j{)}^2}}$$

其中，$X_{ij}$表示第i个样本的第j个特征，${\bar{X}}_j$表示第j个特征的平均值，${\tilde{X}}_{ij}$表示标准化后的特征。

最后，我们可以使用梯度下降或闭式解等方法求解岭回归的参数。

Lasso 回归 (Least Absolute Shrinkage and Selection Operator Regression)

Lasso回归的目标函数如下：

$$\underset{\beta }{\min }||y-X\beta |{|}^2+\lambda ||\beta |{|}_1$$

这里，$X$ 是特征矩阵，$y$ 是目标向量，$\beta$ 是待估参数向量，$\lambda$ 是一个非负的调节参数。

注意到Lasso回归的正则项是$\beta$的L1范数，这使得估计的系数可能会等于零，即某些特征的权重可能被完全剔除，这也是Lasso回归进行特征选择的关键。因为在Lasso回归的目标函数中，$\lambda ||\beta |{|}_1$是一个以原点为顶点的菱形，这就导致优化问题的解可能出现在菱形的角点上。而在Ridge回归的目标函数中，$\lambda ||\beta |{|}^2$是一个以原点为中心的圆形或者超球体，在多维空间中，这个形状没有尖角，优化问题的解不太可能出现在边界上，这就减小了Ridge回归将系数压缩到零的可能性，更可能将系数压缩到一个非常小但非零的值。

Hence, much like best subset selection, the lasso performs feature selection (in an embedded manner). We say that the Lasso yields sparse models, i.e., models that involve only a subset of variables.

因为Lasso回归可能会使一些特征的系数等于零，所以它在处理具有多重共线性问题的数据集时尤为有用。然而，Lasso回归并不能很好地处理特征数量大于样本数量的情况（p > n的情况），在这种情况下，最多只有n个特征能被选择。

Lasso回归也需要进行特征标准化以确保所有特征都在同一尺度上。

Dimension Reduction

主成分分析 (Principal Components Analysis, PCA)

PCA的目标是在高维数据中找到一组新的主成分，它们能最大程度地解释数据的方差。这些新的主成分是原始数据特征的线性组合，原始数据空间的每个方向被替换为一个新的方向，这个新的方向就是主成分。

开始时，我们有一个数据矩阵$X=[X_1,X_2,....,X_p]$，假设它的均值为零。第一个主成分$Z_1$是一个特征的线性组合，它的方差是最大的，即它能最大程度地解释数据的变化。${\varphi }_1$是使得这个线性组合的方差最大的系数向量，我们称其为第一主成分的loadings。

$$Z_1={\varphi }_{11}{X}_1+{\varphi }_{21}{X}_2+\dots +{\varphi }_{p1}{X}_p=\sum _{i=1}^p{\varphi }_{i1}{X}_i=X{\varphi }_1$$

其中 ${\varphi }_1=({\varphi }_{11},\dots ,{\varphi }_{p1}{)}^T$是使得这个线性组合的方差最大的系数向量，我们称其为第一主成分的loadings。正则化意味着loadings的平方和必须为1，即 ${\varphi }_1^T{\varphi }_1=1$。

我们注意到，$Z_1$ 实际上是一个新的向量，其每一个元素都是原始数据的一个观测值在第一主成分上的投影。因此，$Z_1$ 的第 $i$ 个元素可以写为：

$$Z_{i1}=\sum _{j=1}^p{\varphi }_{j1}{x}_{ij}$$

这里，$x_{ij}$ 是数据矩阵 $X$ 的第 $i$ 行第 $j$ 列的元素。

因此，$Z_1$ 的样本方差可以写为：

$$Var(Z_1)=\frac{1}{n-1}\sum _{i=1}^n{Z}_{i1}^2=\frac{1}{n-1}\sum _{i=1}^n(\sum _{j=1}^p{\varphi }_{j1}{x}_{ij}{)}^2$$

在 PCA 中，我们的目标是找到一组加载项 ${\varphi }_{j1}$，使得新向量 $Z_1$ 的样本方差最大。

接下来的问题是，如何找到这个${\varphi }_1=[{\varphi }_{11},{\varphi }_{21},\dots ,{\varphi }_{p1}{]}^T$？实际上，我们需要解决一个带约束的优化问题：我们要找到一个${\varphi }_1$，使得线性组合$Z_1$的方差最大，同时要满足${\varphi }_1$的元素的平方和为1（这是一个正则化的过程，防止${\varphi }_1$取非常大的值）：

$$\underset{{\varphi }_{11},\dots ,{\varphi }_{p1}}{\max }\frac{1}{n-1}\sum _{i=1}^n{\left(\sum _{j=1}^p{\varphi }_{j1}{x}_{ij}\right)}^2\text{ subject to }\sum _{j=1}^p{\varphi }_{j1}^2=1$$

$$\underset{{\varphi }_1}{\max }\frac{1}{n-1}\sum _{i=1}^n(X_i^T{\varphi }_1{)}^2\text{ subject to }{\varphi }_1^T{\varphi }_1=1$$

其中，$X_i$ 是 $X$ 的第 $i$ 行。

由于 $X$ 的每一列的均值都为零，因此 $X$ 的样本协方差矩阵是 $\Sigma =\frac{1}{n-1}{X}^{T}X$。

所以我们有：

$${\varphi }_1^T\Sigma {\varphi }_1={\varphi }_1^T(\frac{1}{n-1}{X}^{T}X){\varphi }_1=\frac{1}{n-1}\sum _{i=1}^n(X_i^T{\varphi }_1{)}^2$$

因此，上述优化问题可以等价于：

$$\underset{{\varphi }_1}{\max }{\varphi }_1^T\Sigma {\varphi }_1\text{ subject to }{\varphi }_1^T{\varphi }_1=1$$

因此，我们可以看到，原始的优化目标实际上等价于最大化 ${\varphi }_1^T\Sigma {\varphi }_1$。所以，通过求解协方差矩阵 $\Sigma$ 的最大特征值对应的特征向量，我们就可以得到使得 $Z_1$ 的样本方差最大化的加载项 ${\varphi }_1$。

Lagrangian函数为：

$$L(\varphi ,\lambda )={\varphi }_1^T\Sigma {\varphi }_1+\lambda (1-{\varphi }_1^T{\varphi }_1)$$

计算偏导数并求解得到：

$$\Sigma {\varphi }_1={\lambda }_1{\varphi }_1$$

$${\varphi }_1^T{\varphi }_1=1$$

这就是特征值方程。得到的特征向量就是loadings。根据上述结果，优化问题变为：

$$\underset{{\varphi }_1,\lambda }{\max }\lambda \text{ subject to }{\varphi }_1^T{\varphi }_1=1,\Sigma {\varphi }_1={\lambda }_1{\varphi }_1$$

因此，$\lambda$ 是 $\Sigma$ 的最大特征值，而 ${\varphi }_1$ 是对应的特征向量。

求解出第一主成分后，我们可以用同样的方法求第二主成分，然后是第三主成分，以此类推。但是，求解第二主成分时，我们除了上述的约束外，还要额外加一个约束：第二主成分要与第一主成分正交，即第二主成分的加载项${\varphi }_2$要满足${\varphi }_2^T{\varphi }_1=0$。

$$\underset{{\varphi }_{12},\dots ,{\varphi }_{p2}}{\max }\frac{1}{n-1}\sum _{i=1}^n{\left(\sum _{j=1}^p{\varphi }_{j2}{x}_{ij}\right)}^2\text{ such that }\sum _{j=1}^p{\varphi }_{j2}^2=1\text{ and }\sum _{j=1}^p{\varphi }_{j1}{\varphi }_{j2}=0$$

或者使用向量表示法：

$$\underset{{\varphi }_2}{\max }{\varphi }_2^T\Sigma {\varphi }_2\text{ such that }{\varphi }_2^T{\varphi }_2=1\text{ and }{\varphi }_2^T{\varphi }_1=0$$

同样，求第三主成分时，除了上述的约束外，还要求第三主成分与第一、第二主成分都正交。其他的主成分以此类推。

最后，我们得到的主成分就形成了一个新的坐标系，这个坐标系下的数据方差最大，同时各个主成分（新的坐标轴）之间是正交的。这样就完成了PCA的过程。

奇异值分解 (SVD)

任何 $n\times m$ 的矩阵 $X$ （假设 $n\ge m$）都可以写成三个矩阵的乘积：

$$X=U\Sigma V^T$$

其中，$U$ 是一个 $n\times m$ 的正交矩阵，$V^T$ 是一个 $m\times m$ 的正交矩阵的转置，$\Sigma$ 是一个对角矩阵，其对角线上的元素是奇异值（非负实数），且按照大小降序排列。

矩阵 $V$ 定义了新的坐标轴，也就是主成分。它的列向量就是主成分载荷。

矩阵 $U$ 是原始数据 $X$ 在新坐标系下的表达，也就是主成分得分。它的第 $i$ 行就是 $X$ 的第 $i$ 行在新坐标系下的坐标。
我们通常只保留前 $k$ 个主成分，也就是只取矩阵 $U$ 和 $V$ 的前 $k$ 列，以及对角矩阵 $\Sigma$ 的前 $k$ 个奇异值。这样，我们就可以将原始数据 $X$ 近似表示为：

$$X\approx U_k{\Sigma }_k{V}_k^T$$

其中，$U_k$ 是 $U$ 的前 $k$ 列，${\Sigma }_k$ 是 $\Sigma$ 的前 $k$ 个奇异值，$V_k^T$ 是 $V^T$ 的前 $k$ 列。

这种方法将数据的维度从 $m$ 降低到了 $k$，大大减少了数据的规模。并且，由于选择的是方差最大的 $k$ 个主成分，所以丢失的信息相对较少，能够对原始数据有一个较好的近似。

例如，考虑图像压缩问题。一个未压缩的灰度图像有 $n\times m$ 个像素，我们需要存储 $n\times m$ 个整数。如果我们使用前 $k$ 个主成分来压缩图像，我们只需要存储：

• 对角矩阵 ${\Sigma }_k$ 中的 $k$ 个奇异值
• 矩阵 $U_k$ 的前 $k$ 列（$n$ 行 $k$ 列）
• 矩阵 $V_k$ 的前 $k$ 列（$m$ 行 $k$ 列）

所以，总共需要存储的数目是 $k(1+n+m)$，压缩比为$\frac{k(1+n+m)}{nm}$。

Scaling

在PCA中，经常会对变量进行中心化处理，也就是移除均值，使每个变量的平均值为零。这是因为PCA是基于方差的，如果我们不移动数据点到均值为0，那么PCA将会更加偏向于那些原始均值更大的变量。

另外，我们也可以通过标准化数据来使所有变量具有标准偏差为1。标准化的目的是为了消除不同变量之间可能存在的量级差异。如果变量的单位不同，其方差会非常不同。因为载荷会更偏向于方差较大的变量，所以如果不进行标准化，那些原本方差较大或者测量单位相对较大的变量将会对主成分分析的结果产生更大的影响。这可能不是我们想要的，因为这可能会使我们忽视那些在原始单位上较小但是在实际影响上较大的变量。

但是，如果所有变量共享相同的单位，那么标准化可能就不必要了，因为没有量纲问题。

性质

• Unique and global solution:

  PCA的结果是确定的，不依赖于初始化或随机种子。

• Ordered components:

  PCA的主成分按照解释方差的多少进行排序，第一主成分解释最多的方差，第二主成分解释次多的方差，依此类推。

• Best low rank approximation to the data:

  数据的最佳低秩近似是指寻找一个秩为 $r$ 的矩阵 $\hat{X}$ 来近似原始矩阵 $X$，使得它们之间的 Frobenius 范数最小。这可以被看作是一个优化问题，即${\min }_{\hat{X}}||X-\hat{X}|{|}_F^2\text{ such that }\text{rank}(\hat{X})\le r$。

  在这个表达式中，$||\cdot |{|}_F$ 是 Frobenius 范数，它定义为矩阵中所有元素平方和的平方根。 $r$ 是我们期望的近似矩阵的秩，它通常远小于 $n$ 和 $p$ （$X$ 的行数和列数）。

  如果我们将数据矩阵 $X$ 投影到由前 $r$ 个主成分载荷向量（即 ${\varphi }_1,\dots ,{\varphi }_r$）定义的子空间，得到的投影矩阵就是最佳的低秩近似。

• Is not the best for non-linear relationships.

t-SNE

t-SNE是一种用于高维数据的可视化的非线性降维方法。与PCA等线性方法不同，t-SNE可以捕捉到更复杂的非线性关系。然而，t-SNE的计算成本较高，且无法像PCA那样提供可以进行推断的主成分。

t-SNE的具体步骤如下：

1. 构造高维对象对的概率分布：t-SNE首先在高维空间中构造一个概率分布，使得相似的对象有较大的概率被选择，而不相似的对象被选择的概率非常小。这一概率分布通常定义为高斯分布，即对于每个点 $i$ ，定义其与点 $j$ 的条件概率 $p_{j|i}$ 如下：

   $$p_{j|i}=\frac{\exp (-||x_i-x_j|{|}^2/2{\sigma }_i^2)}{{\sum }_{k\ne i}\exp (-||x_i-x_k|{|}^2/2{\sigma }_i^2)}$$

   其中，$||x_i-x_j|{|}^2$ 是两点的欧氏距离的平方，${\sigma }_i^2$ 是高斯分布的方差。然后，将这个条件概率对称化，得到点 $i$ 和 $j$ 被同时选择的概率 $p_{ij}$：

   $$p_{ij}=\frac{p_{j|i}+p_{i|j}}{2n}$$

2. 在低维空间中定义相似的概率分布：t-SNE在低维空间中也定义一个类似的概率分布。然而，为了解决在低维空间中的“拥挤问题”，这里采用的是 t 分布而不是高斯分布。假设 $y_i$ 和 $y_j$ 是在低维空间中的两点，定义它们被同时选择的概率 $q_{ij}$ 如下：

   $$q_{ij}=\frac{(1+||y_i-y_j|{|}^2{)}^{-1}}{{\sum }_{k\ne l}(1+||y_k-y_l|{|}^2{)}^{-1}}$$

3. 最小化概率分布之间的 Kullback-Leibler 散度：最后，t-SNE的目标是找到低维空间的点的位置，使得 $q_{ij}$ 尽可能接近 $p_{ij}$。这一目标可以通过最小化 $p_{ij}$ 和 $q_{ij}$ 之间的 Kullback-Leibler 散度来实现：

   $$KL(P||Q)=\sum _{i\ne j}{p}_{ij}\log \frac{p_{ij}}{q_{ij}}$$

   这个最优化问题通常使用梯度下降方法进行求解。

Multidimensional Scaling (MDS)

MDS一种用于在较低维度空间（通常为2维或3维空间）中可视化地表示对象之间近似度（相似性或距离）的技术。与PCA和t-SNE这样的方法不同，MDS只需要距离或非相似性矩阵，而不需要完整的数据集。

MDS通过优化一个损失函数来保留原始数据中的距离，使得在低维空间的投影能尽可能地保留原始数据的距离。已知一个描述数据点之间距离或相似性的矩阵 $D=[d_{ij}]$，寻找一组在低维空间（例如二维或三维）中的投影 $z_1,...,z_k$ ，使得以下的损失函数（也被称为stress function）最小：
在k维空间中，我们有n个点的投影$z_1$, …, $z_n$。我们希望这些投影的欧氏距离尽可能地接近原始的距离
$$S(D,Z)=\sqrt{\frac{{\sum }_{i<j}(d_{ij}-||z_i-z_j||{)}^2}{{\sum }_{i<j}{d}_{ij}^2}}$$

$$S_{LS}(z1,z2,...,zk)=\sqrt{\sum _{i\ne j}(d_{ij}-||z_i-z_j||{)}^2}$$

其中，$d_{ij}$ 是矩阵$D$中的元素，表示数据点 $i$ 和 $j$ 在原始空间中的距离或相似性，$||z_i-z_j||$ 是数据点 $i$ 和 $j$ 在低维空间中的欧氏距离。

解释MDS图的时候，旋转是可以的（因为轴和方向是有一定任意性的），只有相对位置是重要的。我们通常寻找在MDS图中接近的对象，这些对象在原始高维空间中也是相似的。然而，MDS的一个主要挑战是选择合适的低维空间的维度 $k$，这需要根据问题的具体情况来确定，或者通过一些启发式的方法（如Scree图中的肘部法则）来估计。