有时数据集会变得非常大，甚至可能非常高维。然而，在许多情况下，数据本身是稀疏的——也就是说，尽管有许多特征，但任何给定的样本只有少数非零特征被观察到。在这种情况下，通过稀疏矩阵数据结构可以更有效地表示数据，从而节省内存使用。找到直接适用于这种稀疏数据的降维技术可能很困难——通常人们会应用一种基本的线性技术，例如来自sklearn的TruncatedSVD（它确实接受稀疏矩阵输入），以将数据转换为适合其他更高级降维技术的格式。在UMAP的情况下，这不是必要的——UMAP可以直接在稀疏矩阵输入上运行。本教程将通过几个示例来展示如何做到这一点。首先，我们需要加载一些库。显然我们需要numpy，但我们还将使用scipy.sparse，它提供了稀疏矩阵数据结构。我们的一个示例将是纯数学的，我们将为此使用sympy；另一个示例是基于文本的，我们将为此使用sklearn（特别是sklearn.feature_extraction.text）。除此之外，我们还需要umap和绘图工具。

import numpy as np
import scipy.sparse
import sympy
import sklearn.datasets
import sklearn.feature_extraction.text
import umap
import umap.plot
import matplotlib.pyplot as plt
%matplotlib inline

一个数学示例

要开始，我们需要一个所有素数的列表。幸运的是，我们有sympy可以使用，并且我们可以通过一次调用primerange快速获取该信息。我们还需要一个字典，将不同的素数映射到它们在我们数据结构中对应的列号；实际上，我们只是对素数进行枚举。

primes = list(sympy.primerange(2, 110000))
prime_to_column = {p:i for i, p in enumerate(primes)}

现在我们需要以一种可以轻松放入稀疏矩阵的格式来构建我们的数据。在这一点上，了解一些关于稀疏矩阵数据结构的背景知识是有用的。为此，我们将使用所谓的“LIL”格式。 LIL是“列表的列表”的缩写，因为数据在内部就是这样存储的。有一个包含所有行的列表，每一行都存储为一个列表，给出非零项的列索引。为了存储数据值，有一个并行的结构，包含与给定行和列对应的项的值。

与此同时，我们将构建相应的值结构以插入矩阵中。由于我们只关心可除性，这将在每个非零条目中简单地表示为1，因此我们可以为每一行添加一个适当长度的1的列表。

%%time
lil_matrix_rows = []
lil_matrix_data = []
for n in range(100000):
    prime_factors = sympy.primefactors(n)
    lil_matrix_rows.append([prime_to_column[p] for p in prime_factors])
    lil_matrix_data.append([1] * len(prime_factors))

CPU times: user 2.07 s, sys: 26.4 ms, total: 2.1 s
Wall time: 2.1 s

现在我们需要将其转换为稀疏矩阵。幸运的是，scipy.sparse 包使这变得容易，而且我们已经以一种相当有用的结构构建了数据。首先，我们创建一个正确格式（LIL）和正确形状（与我们生成的行数相同，列数与素数数量相同）的稀疏矩阵。这本质上只是一个空矩阵。我们可以通过将 rows 属性设置为我们生成的行，并将 data 属性设置为相应的值结构（全部为1）来修复这个问题。结果是一个稀疏矩阵数据结构，可以轻松操作并转换为其他稀疏矩阵格式。

factor_matrix = scipy.sparse.lil_matrix((len(lil_matrix_rows), len(primes)), dtype=np.float32)
factor_matrix.rows = np.array(lil_matrix_rows)
factor_matrix.data = np.array(lil_matrix_data)
factor_matrix

<100000x10453 sparse matrix of type '<class 'numpy.float32'>'
    with 266398 stored elements in LInked List format>

如你所见，我们有一个包含100000行和超过10000列的矩阵。 如果我们将其存储为numpy数组，它将占用大量内存。然而，实际上只有大约260000个条目是非零的，而这正是我们真正需要存储的，使其更加紧凑。

%%time
mapper = umap.UMAP(metric='cosine', random_state=42, low_memory=True).fit(factor_matrix)

CPU times: user 9min 36s, sys: 6.76 s, total: 9min 43s
Wall time: 9min 7s

这很简单！但它真的有效吗？我们可以轻松地绘制结果：

umap.plot.points(mapper, values=np.arange(100000), theme='viridis')
_images/sparse_11_1.png

这看起来与John Williamson 得到的结果非常一致，前提是我们只使用了100,000个整数而不是1,000,000个，以确保大多数用户应该能够运行这个示例（完整的百万可能需要一个大型内存计算节点）。所以看起来这个工作得很好。接下来的问题是，我们是否可以使用transform功能将新数据映射到这个空间。为了测试这一点，我们需要更多的数据。幸运的是，还有更多的整数。我们将抓取接下来的10,000个并将它们放入稀疏矩阵中，就像我们对前100,000个所做的那样。

%%time
lil_matrix_rows = []
lil_matrix_data = []
for n in range(100000, 110000):
    prime_factors = sympy.primefactors(n)
    lil_matrix_rows.append([prime_to_column[p] for p in prime_factors])
    lil_matrix_data.append([1] * len(prime_factors))

CPU times: user 214 ms, sys: 1.99 ms, total: 216 ms
Wall time: 222 ms

new_data = scipy.sparse.lil_matrix((len(lil_matrix_rows), len(primes)), dtype=np.float32)
new_data.rows = np.array(lil_matrix_rows)
new_data.data = np.array(lil_matrix_data)
new_data

<10000x10453 sparse matrix of type '<class 'numpy.float32'>'
    with 27592 stored elements in LInked List format>

为了映射我们生成的新数据，我们可以简单地将其传递给训练模型的transform方法。这有点慢，但确实有效。

new_data_embedding = mapper.transform(new_data)

我们可以绘制结果。由于这次我们只得到了点的位置（而不是模型），我们将不得不使用matplotlib进行绘图。

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111)
plt.scatter(new_data_embedding[:, 0], new_data_embedding[:, 1], s=0.1, c=np.arange(10000), cmap='viridis')
ax.set(xticks=[], yticks=[], facecolor='black');
_images/sparse_18_0.png

在这种情况下，颜色比例不同，但你可以看到数据已经被映射到与原始嵌入中看到的各种结构相对应的位置。因此，即使对于大型稀疏数据，我们也可以嵌入数据，甚至可以将新数据添加到嵌入中。

一个文本分析示例

在这个例子中，我们将使用经典的20newsgroups数据集，这是从旧的NNTP新闻组系统中抽取的新闻组消息样本，涵盖了20个不同的新闻组。sklearn.datasets模块可以轻松获取这些数据，实际上，我们可以获取一个预向量化的版本来省去我们自己运行CountVectorizer的麻烦。我们将获取训练集和测试集以供后续使用。

news_train = sklearn.datasets.fetch_20newsgroups_vectorized(subset='train')
news_test = sklearn.datasets.fetch_20newsgroups_vectorized(subset='test')

如果我们查看我们拉取的实际数据，我们会看到 sklearn 已经运行了一个 CountVectorizer 并以稀疏 矩阵格式生成了数据。

news_train.data

<11314x130107 sparse matrix of type '<class 'numpy.float64'>'
    with 1787565 stored elements in Compressed Sparse Row format>

在这种情况下，稀疏矩阵格式的价值立即显现出来；虽然只有11,000个样本，但有130,000个特征！如果数据存储在标准的numpy数组中，我们将占用10GB的内存！而且大部分内存只是反复存储数字零。在稀疏矩阵格式中，它很容易适应大多数机器的内存。这种数据维度在文本工作负载中非常常见。

然而，原始计数并不理想，因为像“the”和“and”这样的常见词会主导大多数文档的计数，而对文档实际内容的信息贡献却很少。我们可以通过使用sklearn中的TfidfTransformer来纠正这一点，它将数据转换为TF-IDF格式。关于TF-IDF所做的转换有很多种理解方式，但我喜欢直观地将其理解为：一个词的信息量可以（大致）与其频率的负对数成正比；一个词使用得越频繁，它携带的信息量就越少，而不常见的词则携带更多的信息。TF-IDF的作用可以被视为根据与该列相关的词的信息量重新加权列。因此，像“the”和“and”这样的常见词会被降权，因为它们携带的关于文档的信息较少，而不常见的词则被认为更重要，其相关的列会被加权。我们可以将这种转换应用于训练集和测试集（使用在训练集上训练的相同转换器）。

tfidf = sklearn.feature_extraction.text.TfidfTransformer(norm='l1').fit(news_train.data)
train_data = tfidf.transform(news_train.data)
test_data = tfidf.transform(news_test.data)

结果仍然是一个稀疏矩阵，因为TF-IDF根本不会改变零元素，也不会改变特征的数量。

train_data

<11314x130107 sparse matrix of type '<class 'numpy.float64'>'
    with 1787565 stored elements in Compressed Sparse Row format>

现在我们需要将这个非常高维的数据传递给UMAP。与一些其他非线性降维技术不同，我们不需要先应用PCA将数据降到合理的维度；我们也不需要使用其他技术来减少数据以便能够表示为密集的numpy数组；我们可以直接在130,000维的稀疏矩阵上工作。

%%time
mapper = umap.UMAP(metric='hellinger', random_state=42).fit(train_data)

CPU times: user 8min 40s, sys: 3.07 s, total: 8min 44s
Wall time: 8min 43s

umap.plot.points(mapper, labels=news_train.target)
_images/sparse_31_1.png

我们可以看到，即使直接从130,000维空间降到仅2维，UMAP在分离许多不同的新闻组方面做得相当不错。

我们现在可以尝试使用transform方法将测试数据添加到相同的空间中。

test_embedding = mapper.transform(test_data)

fig = plt.figure(figsize=(12,12))
ax = fig.add_subplot(111)
plt.scatter(test_embedding[:, 0], test_embedding[:, 1], s=1, c=news_test.target, cmap='Spectral')
ax.set(xticks=[], yticks=[]);
_images/sparse_35_0.png