R 中的数据挖掘算法/降维/奇异值分解

在本章中，我们将探讨奇异值分解（SVD），这是一种矩阵分解方法，利用线性代数知识来进行这样的分解。

奇异值分解

在数据挖掘中，此算法可用于通过显示重要维度的数量来更好地理解数据库，并通过减少数据挖掘过程中使用的属性数量来简化数据库。这种减少消除了从线性代数角度来看线性相关的非必要数据。例如，假设一个数据库包含一个字段，该字段存储多个样本的水温，另一个字段存储其状态（固体、液体或气体）。很容易看出，第二个字段依赖于第一个字段，因此，SVD 可以轻松地向我们表明它对于分析并不重要。

主成分分析 (PCA) 是 SVD 的一个特例。

算法

SVD 是对具有$n$ 行和 $d$ 列的实数或复数矩阵 $X$ 的因式分解。

$X=UAV^{\top }$

其中 $U$ 是一个维度为 $n\times n$ 的矩阵，$V$ 是另一个维度为 $d\times d$ 的矩阵，$A$ 是一个维度为 $n\times d$ 的矩阵，与 $X$ 的维度相同。

此外

$U^{\top }U=I_{n}$ 和 $V^{\top }V=I_{d}$

其中 $I_{n}$ 和 $I_{d}$ 是单位矩阵，大小分别为 $n$ 和 $d$.

矩阵 $U$ 的列是矩阵 $X$ 的 **左奇异向量**，矩阵 $V$ （或 $V^{\top }$ 的行）的列是 **右奇异向量**。

矩阵 $A$ 是对角矩阵，其对角线上的值是矩阵 $X$ 的 **奇异值**。矩阵 $A$ 中每行上的奇异值永远不会小于其下方行的值。所有奇异值都大于 $0$.

计算 SVD 就是找到 $XX^{\top }$ 和 $X^{\top }X$ 的特征值和特征向量。 $X^{\top }X$ 的特征向量是 $V$ 的列，而 $XX^{\top }$ 的特征向量是 $U$ 的列。矩阵 $X$ 的奇异值（在矩阵 $A$ 的对角线上）是 $XX^{\top }$ 和 $X^{\top }X$ 共有的正特征值的平方根。

如果 $XX^{\top }$ 和 $X^{\top }X$ 具有相同数量的特征值，则 $X$ 是一个方阵；否则，特征值较少的矩阵的特征值是特征值较多的矩阵的特征值。因此，$X$ 的奇异值是矩阵的特征值，位于 $XX^{\top }$ 和 $X^{\top }X$ 中，特征值较少的矩阵中。

矩阵的奇异值数量等于该矩阵的秩，即矩阵中线性无关的列或行的数量。秩不大于 $\min(n,d)$，因为这是矩阵对角线上元素的数量。奇异值是矩阵 $A$ 对角线上的元素。正奇异值的数量等于矩阵的秩。

因此，算法如下

1. 计算 $XX^{\top }$ 通常。
2. 计算 $XX^{\top }$ 的特征值和特征向量通常。
3. 计算 $X^{\top }X$ 。
4. 计算 $X^{\top }X$ 的特征值和特征向量。
5. 计算 $XX^{\top }$ 和 $X^{\top }X$ 的公共正特征值的平方根。
6. 最后，将计算出的值分配给 $U$，$V$ 和 $A$ 。

示例

X = ${\begin{bmatrix}1&1&0&2\\1&1&2&0\\2&0&1&1\end{bmatrix}}$

XX^T = ${\begin{bmatrix}6&2&4\\2&6&4\\4&4&6\end{bmatrix}}$

XX^T 的特征值为

12.744563, 4.000000 和 1.255437

XX^T 的特征向量为

1. 0.5417743, 0.5417743, 0.6426206
2. 0.7071068, -0.7071068, -2.144010 x 10^-16
3. 0.4544013, 0.4544013, -0.7661846

X^TX = ${\begin{bmatrix}6&2&4&4\\2&2&2&2\\4&2&5&1\\4&2&1&5\end{bmatrix}}$

X^TX 的特征值为

12.744563, 4.000000, 1.255437 和 (5.940557 x 10^-18)

X^TX 的特征向量为

1. 0.6635353, 0.3035190, 0.4835272, 0.4835272
2. 0.0000000, 0.5181041 x 10^-16, -0.7071068, 0.7071068
3. -0.5565257, 0.8110957, 0.1272850, 0.1272850
4. 0.5000000, 0.5000000, -0.5000000, -0.5000000

X 的奇异值为 X^TX 和 XX^T 的共同特征值的平方根

${\sqrt {12.744563}}=3.569953$

${\sqrt {4.000000}}=2.000000$

${\sqrt {1.255437}}=1.120463$

因此

A = ${\begin{bmatrix}3.569953&0&0&0\\0&2.000000&0&0\\0&0&1.120463&0\end{bmatrix}}$

最后，X 被分解为

X = UAV^T = ${\begin{bmatrix}0.5417743&0.7071068&0.4544013\\0.5417743&-0.7071068&0.4544013\\0.6426206&-2.144010\cdot 10^{-16}&-0.7661846\end{bmatrix}}{\begin{bmatrix}3.569953&0&0&0\\0&2.000000&0&0\\0&0&1.120463&0\end{bmatrix}}{\begin{bmatrix}0.6635353&0.000000&-0.5565257&0.5\\0.3035190&5.181041\cdot 10^{-16}&0.8110957&0.5\\0.4835272&-0.7071068&0.1272850&-0.5\\0.4835272&0.7071068&0.1272850&-0.5\end{bmatrix}}$

实现

R 内置了一个名为'svd()'的函数用于计算 SVD。

默认情况下，它接收 R 的 原生矩阵 作为参数，并返回一个包含 U、A 和 V 的 数据框。

参数

如果不需要计算所有奇异值和向量，可以告诉 svd() 需要多少个元素。

这可以通过将这些值分别分配给 **nu** 和 **nv** 来实现，它们分别代表所需的左奇异向量和右奇异向量数量。

例如，为了只计算一半的向量，可以执行以下操作

svd(X, nu = min(nrow(X), ncol(X)) / 2, nv = min(nrow(X), ncol(X))/2)

返回的对象

s = svd(X)

考虑到

X = UAV^T

返回的对象是一个包含三个字段的数据结构

s$d 是包含 X 的奇异值的向量，它来自矩阵 A 的对角线。

s$u 是一个矩阵，它的列包含 X 的左奇异向量。它的行数与 X 的行数相同，它的列数与传递给参数 nu 的数字相同。注意，如果 nu 为 0，则不会创建此矩阵。

s$v 是一个矩阵，它的列包含 X 的右奇异向量。它的行数与 X 的列数相同，它的列数与传递给参数 nv 的数字相同。注意，如果 nv 为 0，则不会创建此矩阵。

例子

在 R 终端或 R 程序中执行以下命令序列，并查看结果

示例 1

dat <- seq(1,240,2)

X <- matrix(dat,ncol=12)

s <- svd(X)

A <- diag(s$d)

s$u %*% A %*% t(s$v) #  X = U A V'

示例 2

dat <- seq(1,240,2)

X <- matrix(dat,ncol=12)

s <- svd(X, nu = nrow(X), nv = ncol(X))

A <- diag(s$d)

A <- cbind(A, 0)  # Add two columns with zero, in order to A have the same dimensions of X.

A <- cbind(A, 0)

s$u %*% A %*% t(s$v) #  X = U A V'

案例研究

场景

为了更好地可视化 SVD 的结果，在这个案例研究中，我们将使用一个表示图像的矩阵，这样任何对矩阵的更改都可以轻松观察到。

输入数据

要使用 R 处理图像，应安装 rimage 包

> install.packages("rimage")

然后将图像加载到 R 中，将其转换为灰度图像。这样，我们将得到一个二进制矩阵，其中 1 表示白色，0 表示黑色。

library(rimage)
tux <- read.jpeg('tux.jpeg')
tux <- imagematrix(tux,type='grey')

我们可以看到此导入的结果

plot(tux)

为了帮助我们进行降维，让我们创建一个小的辅助函数，该函数将接收我们的 tux 和我们想要的维度数，并返回我们的新的 tux。

reduce <- function(A,dim) {
    #Calculates the SVD
    sing <- svd(A)

    #Approximate each result of SVD with the given dimension
    u<-as.matrix(sing$u[, 1:dim])
    v<-as.matrix(sing$v[, 1:dim])
    d<-as.matrix(diag(sing$d)[1:dim, 1:dim])

    #Create the new approximated matrix
    return(imagematrix(u%*%d%*%t(v),type='grey'))
}

执行和输出

现在我们有了矩阵，让我们看看它有多少奇异值

tux_d <- svd(tux)
length(tux_d$d)
[1] 335

因此，这个矩阵有 335 个奇异值。让我们首先尝试将其减少到只有一个奇异值

plot(reduce(tux,1))

正如我们所看到的，这种近似值从我们的矩阵中删除了许多有用的信息。如果它是数据库，我们肯定已经丢失了重要的数据。

让我们尝试使用 10% (35) 个奇异值

plot(reduce(tux,35))

分析

仅使用 10% 的真实数据，我们就能创建出真实数据的非常好的近似值。此外，通过这种方法，我们只需使用最重要的奇异值就可以去除噪声和线性相关元素。这在数据挖掘中非常有用，因为很难确定数据库是否“干净”，如果不是，哪些元素对我们的分析有用或无用。

参考文献

• 数据挖掘算法基础 - Mohammed J. Zaki，Wagner Meira Jr. - 第 7.4 章
• 奇异值分解
• http://www.ats.ucla.edu/stat/r/code/svd_demos.htm
• http://stat.ethz.ch/R-manual/R-patched/library/base/html/svd.html