我有一个大的、稀疏的二进制矩阵(大约39,000 x 14,000; 大多数行只有一个"1"条目)。我想把相似的行聚类在一起，但我的初始计划需要太长时间才能完成:
d <- dist(inputMatrix, method="binary")
hc <- hclust(d, method="complete")
第一步尚未完成，因此我不确定第二步会如何处理。有哪些方法可以在R中高效地对大型、稀疏、二元矩阵中的相似行进行分组？
在R中对一个大的、非常稀疏的二元矩阵进行聚类
9
- Matt LaFave
6
3个回答
20
我写了一些Rcpp代码和R代码，可以比
dist(x,method="binary")函数快约80倍地计算二进制矩阵的二进制/Jaccard距离。它将输入矩阵转换为原始矩阵，原始矩阵是输入矩阵的转置（以便位模式在内部正确排序）。然后，使用一些C++代码处理数据作为64位无符号整数以提高速度。两个向量x和y的Jaccard距离等于x ^ y / (x | y)，其中^是异或运算符。如果xor或or的结果非零，则使用汉明重量计算来计算设置的位数。
我已经在github上放置了代码：https://github.com/NikNakk/binaryDist/并在下面复制了两个文件。我确认对于一些随机数据集，结果与dist(x,method="binary")相同。
在一个由39000行和14000列组成，每行1-5个元素的数据集中，它大约需要11分钟。输出距离矩阵为5.7 GB。
bDist.cpp
#include <Rcpp.h>
using namespace Rcpp;
//countBits function taken from https://en.wikipedia.org/wiki/Hamming_weight#Efficient_implementation
const uint64_t m1 = 0x5555555555555555; //binary: 0101...
const uint64_t m2 = 0x3333333333333333; //binary: 00110011..
const uint64_t m4 = 0x0f0f0f0f0f0f0f0f; //binary: 4 zeros, 4 ones ...
const uint64_t h01 = 0x0101010101010101; //the sum of 256 to the power of 0,1,2,3...
int countBits(uint64_t x) {
x -= (x >> 1) & m1; //put count of each 2 bits into those 2 bits
x = (x & m2) + ((x >> 2) & m2); //put count of each 4 bits into those 4 bits
x = (x + (x >> 4)) & m4; //put count of each 8 bits into those 8 bits
return (x * h01)>>56; //returns left 8 bits of x + (x<<8) + (x<<16) + (x<<24) + ...
}
// [[Rcpp::export]]
int countBitsFromRaw(RawVector rv) {
uint64_t* x = (uint64_t*)RAW(rv);
return(countBits(*x));
}
// [[Rcpp::export]]
NumericVector bDist(RawMatrix mat) {
int nr(mat.nrow()), nc(mat.ncol());
int nw = nr / 8;
NumericVector res(nc * (nc - 1) / 2);
// Access the raw data as unsigned 64 bit integers
uint64_t* data = (uint64_t*)RAW(mat);
uint64_t a(0);
// Work through each possible combination of columns (rows in the original integer matrix)
for (int i = 0; i < nc - 1; i++) {
for (int j = i + 1; j < nc; j++) {
uint64_t sx = 0;
uint64_t so = 0;
// Work through each 64 bit integer and calculate the sum of (x ^ y) and (x | y)
for (int k = 0; k < nw; k++) {
uint64_t o = data[nw * i + k] | data[nw * j + k];
// If (x | y == 0) then (x ^ y) will also be 0
if (o) {
// Use Hamming weight method to calculate number of set bits
so = so + countBits(o);
uint64_t x = data[nw * i + k] ^ data[nw * j + k];
if (x) {
sx = sx + countBits(x);
}
}
}
res(a++) = (double)sx / so;
}
}
return (res);
}
R语言源代码
library("Rcpp")
library("plyr")
sourceCpp("bDist.cpp")
# Converts a binary integer vector into a packed raw vector,
# padding out at the end to make the input length a multiple of packWidth
packRow <- function(row, packWidth = 64L) {
packBits(as.raw(c(row, rep(0, (packWidth - length(row)) %% packWidth))))
}
as.PackedMatrix <- function(x, packWidth = 64L) {
UseMethod("as.PackedMatrix")
}
# Converts a binary integer matrix into a packed raw matrix
# padding out at the end to make the input length a multiple of packWidth
as.PackedMatrix.matrix <- function(x, packWidth = 64L) {
stopifnot(packWidth %% 8 == 0, class(x) %in% c("matrix", "Matrix"))
storage.mode(x) <- "raw"
if (ncol(x) %% packWidth != 0) {
x <- cbind(x, matrix(0L, nrow = nrow(x), ncol = packWidth - (ncol(x) %% packWidth)))
}
out <- packBits(t(x))
dim(out) <- c(ncol(x) %/% 8, nrow(x))
class(out) <- "PackedMatrix"
out
}
# Converts back to an integer matrix
as.matrix.PackedMatrix <- function(x) {
out <- rawToBits(x)
dim(out) <- c(nrow(x) * 8L, ncol(x))
storage.mode(out) <- "integer"
t(out)
}
# Generates random sparse data for testing the main function
makeRandomData <- function(nObs, nVariables, maxBits, packed = FALSE) {
x <- replicate(nObs, {
y <- integer(nVariables)
y[sample(nVariables, sample(maxBits, 1))] <- 1L
if (packed) {
packRow(y, 64L)
} else {
y
}
})
if (packed) {
class(x) <- "PackedMatrix"
x
} else {
t(x)
}
}
# Reads a binary matrix from file or character vector
# Borrows the first bit of code from read.table
readPackedMatrix <- function(file = NULL, text = NULL, packWidth = 64L) {
if (missing(file) && !missing(text)) {
file <- textConnection(text)
on.exit(close(file))
}
if (is.character(file)) {
file <- file(file, "rt")
on.exit(close(file))
}
if (!inherits(file, "connection"))
stop("'file' must be a character string or connection")
if (!isOpen(file, "rt")) {
open(file, "rt")
on.exit(close(file))
}
lst <- list()
i <- 1
while(length(line <- readLines(file, n = 1)) > 0) {
lst[[i]] <- packRow(as.integer(strsplit(line, "", fixed = TRUE)[[1]]), packWidth = packWidth)
i <- i + 1
}
out <- do.call("cbind", lst)
class(out) <- "PackedMatrix"
out
}
# Wrapper for the C++ code which
binaryDist <- function(x) {
if (class(x) != "PackedMatrix") {
x <- as.PackedMatrix(x)
}
dst <- bDist(x)
attr(dst, "Size") <- ncol(x)
attr(dst, "Diag") <- attr(dst, "Upper") <- FALSE
attr(dst, "method") <- "binary"
attr(dst, "call") <- match.call()
class(dst) <- "dist"
dst
}
x <- makeRandomData(2000, 400, maxBits = 5, packed = TRUE)
system.time(bd <- binaryDist(x))
从原始答案：
其他需要考虑的事情是，在比较两行单个值时进行一些预过滤，因为距离将在重复项时为0，而在任何其他可能性时为1。
几个相对简单的选项，可能更快，不需要太多代码，包括 vegan 包中的 vegdist 函数和 amap 包中的 Dist 函数。如果您有多个核心并利用它支持并行化的优势，后者可能只会更快。
- Nick Kennedy
5
1可以，谢谢！还有一件事 -
bd 与 dist() 的输出略有不同，因为 bd 缺少 "label" 字段。看起来通过 names(cutClusterObject) <- rownames(x) 在聚类后将它们添加回去是有效的，但您认为有办法在 bd 对象中保留标签吗？ - Matt LaFave@MattLaFave 说实话，最简单的方法就是按照你提出的建议去做 - 如果需要的话，这可以添加到上面的代码中。 - Nick Kennedy
这太棒了。谢谢你的发布！ - aaronwolen
1我在我的研究中使用了你的代码，我想在论文中引用它。你能告诉我你希望如何引用吗？ - Shamsa
@Shamsa 谢谢，很高兴听到它对你有用！这段代码除了在这里之外没有发布在其他地方，所以最好引用这个 Stack Overflow 帖子和我的名字。如果您愿意这样做，在您准备与世界分享论文时，我们会非常感兴趣看到它。 - Nick Kennedy
6
这需要很长时间计算的原因是调用
一种更快的方法是k-means聚类，因为它不需要预先计算距离矩阵；在每次迭代中，您只需要进行39000*k次距离计算，其中k是簇的数量。要获得类似于Jaccard指数（如果相同则为0，如果没有索引重合则为1，如果有但不是所有索引重合则在两者之间）的成对距离，您可以通过将每行
dist正在计算并存储超过7.6亿个成对距离。如果您的数据存储得稀疏，这将需要很长时间和大量存储空间。如果您的数据没有存储得稀疏，则每次距离计算至少需要14,000次操作，总操作次数超过1千万亿！一种更快的方法是k-means聚类，因为它不需要预先计算距离矩阵；在每次迭代中，您只需要进行39000*k次距离计算，其中k是簇的数量。要获得类似于Jaccard指数（如果相同则为0，如果没有索引重合则为1，如果有但不是所有索引重合则在两者之间）的成对距离，您可以通过将每行
x除以sqrt(2*sum(x^2))来实现。例如，如果您有以下输入矩阵：(mat <- rbind(c(1, 0, 0, 0, 0), c(0, 0, 0, 1, 1)))
# [,1] [,2] [,3] [,4] [,5]
# [1,] 1 0 0 0 0
# [2,] 0 0 0 1 1
假设矩阵中仅包含二进制值，则规范化版本应为（如果不是这种情况，则使用rowSums（mat ^ 2））：
(mat.norm <- mat / sqrt(2*rowSums(mat)))
# [,1] [,2] [,3] [,4] [,5]
# [1,] 0.7071068 0 0 0.0 0.0
# [2,] 0.0000000 0 0 0.5 0.5
这两个观测值（它们没有共同的指标）之间的欧几里得距离为1，恰好对应于此案例的Jaccard距离。
dist(mat.norm, "euclidean")
# 1
# 2 1
此外，相同的观测明显具有欧几里得距离为0，这也对应于Jaccard距离。
- josliber
2
您是否有重复行？不需要计算它们的距离两次。
所有只有一个1的行都与所有在不同位置上只有一个1的行完全不同。
因此，在这样的数据上运行聚类没有意义。输出结果是可预测的，并且归结为找到1。
尝试仅限于具有多个1的对象的数据集。除非您能够在这些数据上获得有趣的结果，否则无需继续进行。二进制数据信息太少。
- Has QUIT--Anony-Mousse
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接
method设置为“binary”时，这就是dist使用的距离。 - Matt LaFave