1

我有一个非常大的数据集,所以我试图用下面的一个小例子来总结我的问题。

假设我有一个名为 X 的 3X3 矩阵,列名为 a、b 和 c。

X = (1, 10, 0.1,
     2, 20, 0.2,
     3, 30, 0.3)

wherea = c(1, 2, 3)给出要重复的次数,b = c(10, 20, 30)给出要重复的实际值,如果 in 的次数小于 4(矩阵 Y 的列数)c = c(0.1, 0.2, 0.3),则给出要填写的值。a

我的目标是生成一个 3X4 矩阵 Y,应该是这样的

Y = (10, 0.1, 0.1, 0.1,
     20,  20, 0.2, 0.2,
     30,  30,  30, 0.3)

我知道可能有很多方法可以做这个例子,但由于我的真实数据非常大(X 有 100 万行,Y 有 480 列),我真的必须在没有循环的情况下这样做(比如 480 次迭代)。我已经尝试使用该功能rep,但仍然无法做到这一点。

4

2 回答 2

4

输出矩阵的每一行都可以通过对函数的一次调用来计算rep,使整个操作成为 1-liner:

t(apply(X, 1, function(x) rep(x[2:3], c(x[1], 4-x[1]))))
#      [,1] [,2] [,3] [,4]
# [1,]   10  0.1  0.1  0.1
# [2,]   20 20.0  0.2  0.2
# [3,]   30 30.0 30.0  0.3

您说您计划创建一个 1e6 x 480 矩阵,希望它适合您的系统内存。但是,您可能无法在不耗尽系统内存的情况下将其推得更大。

于 2015-05-11T21:56:37.530 回答
1

解决方案

这并不容易,但我想出了一种方法来完成这项任务,使用单个矢量化调用rep(),加上一些脚手架代码:

XR <- 3;
YC <- 4;
X <- matrix(c(1:XR%%(YC+1),seq(10,by=10,length.out=XR),seq(0.1,by=0.1,length.out=XR)),XR,dimnames=list(NULL,c('rep','val','fill')));
X;
##      rep val fill
## [1,]   1  10  0.1
## [2,]   2  20  0.2
## [3,]   3  30  0.3
Y <- matrix(rep(t(X[,c('val','fill')]),times=c(rbind(X[,'rep'],YC-X[,'rep']))),XR,byrow=T);
Y;
##      [,1] [,2] [,3] [,4]
## [1,]   10  0.1  0.1  0.1
## [2,]   20 20.0  0.2  0.2
## [3,]   30 30.0 30.0  0.3

(次要点:我选择将列名分配rep val fillX,而不是a b c问题中指定的,并且我在索引时在我的解决方案中使用了这些列名X(而不是使用数字索引),因为我通常更喜欢最大化人类可读性在可能的情况下,但这个细节对于解决方案的正确性和性能可以忽略不计。)

表现

这实际上比@josilber 的解决方案具有显着的性能优势,因为他使用apply()which 在矩阵的行上进行内部循环(在 R 语言中传统上称为“隐藏循环”),而我的解决方案的核心是对rep(). 我这样说不是为了敲@josilber 的解决方案,这是一个很好的解决方案(我什至给了他一个赞成票!);这不是解决这个问题的最佳解决方案。

这是使用您在问题中指出的大量参数的性能优势演示:

XR <- 1e6;
YC <- 480;
X <- matrix(c(1:XR%%(YC+1),seq(10,by=10,length.out=XR),seq(0.1,by=0.1,length.out=XR)),XR,dimnames=list(NULL,c('rep','val','fill')));
X;
##        rep  val fill
##   [1,]   1   10  0.1
##   [2,]   2   20  0.2
##   [3,]   3   30  0.3
##   [4,]   4   40  0.4
##   [5,]   5   50  0.5
##   [6,]   6   60  0.6
##   [7,]   7   70  0.7
##   [8,]   8   80  0.8
##   [9,]   9   90  0.9
##  [10,]  10  100  1.0
##  [11,]  11  110  1.1
##  [12,]  12  120  1.2
##  [13,]  13  130  1.3
##
## ... (snip) ...
##
## [477,] 477 4770 47.7
## [478,] 478 4780 47.8
## [479,] 479 4790 47.9
## [480,] 480 4800 48.0
## [481,]   0 4810 48.1
## [482,]   1 4820 48.2
## [483,]   2 4830 48.3
## [484,]   3 4840 48.4
## [485,]   4 4850 48.5
## [486,]   5 4860 48.6
## [487,]   6 4870 48.7
## [488,]   7 4880 48.8
## [489,]   8 4890 48.9
## [490,]   9 4900 49.0
## [491,]  10 4910 49.1
## [492,]  11 4920 49.2
##
## ... (snip) ...
##
## [999986,] 468  9999860  99998.6
## [999987,] 469  9999870  99998.7
## [999988,] 470  9999880  99998.8
## [999989,] 471  9999890  99998.9
## [999990,] 472  9999900  99999.0
## [999991,] 473  9999910  99999.1
## [999992,] 474  9999920  99999.2
## [999993,] 475  9999930  99999.3
## [999994,] 476  9999940  99999.4
## [999995,] 477  9999950  99999.5
## [999996,] 478  9999960  99999.6
## [999997,] 479  9999970  99999.7
## [999998,] 480  9999980  99999.8
## [999999,]   0  9999990  99999.9
## [1e+06,]    1 10000000 100000.0
josilber <- function() t(apply(X,1,function(x) rep(x[2:3],c(x[1],YC-x[1]))));
bgoldst <- function() matrix(rep(t(X[,c('val','fill')]),times=c(rbind(X[,'rep'],YC-X[,'rep']))),XR,byrow=T);
system.time({ josilber(); });
##    user  system elapsed
##  65.719   3.828  71.623
system.time({ josilber(); });
##    user  system elapsed
##  60.375   2.609  66.724
system.time({ bgoldst(); });
##    user  system elapsed
##   5.422   0.593   6.033
system.time({ bgoldst(); });
##    user  system elapsed
##   5.203   0.797   6.002

只是为了证明@josilber 和我得到了完全相同的结果,即使对于这么大的输入:

identical(bgoldst(),josilber());
## [1] TRUE

解释

现在我将尝试解释解决方案的工作原理。为了解释,我将使用以下输入:

XR <- 6;
YC <- 4;
X <- matrix(c(1:XR%%(YC+1),seq(10,by=10,length.out=XR),seq(0.1,by=0.1,length.out=XR)),XR,dimnames=list(NULL,c('rep','val','fill')));
X;
##      rep val fill
## [1,]   1  10  0.1
## [2,]   2  20  0.2
## [3,]   3  30  0.3
## [4,]   4  40  0.4
## [5,]   0  50  0.5
## [6,]   1  60  0.6

解决方案是:

Y <- matrix(rep(t(X[,c('val','fill')]),times=c(rbind(X[,'rep'],YC-X[,'rep']))),XR,byrow=T);
Y;
##      [,1] [,2] [,3] [,4]
## [1,] 10.0  0.1  0.1  0.1
## [2,] 20.0 20.0  0.2  0.2
## [3,] 30.0 30.0 30.0  0.3
## [4,] 40.0 40.0 40.0 40.0
## [5,]  0.5  0.5  0.5  0.5
## [6,] 60.0  0.6  0.6  0.6

在高层次上,解决方案是围绕形成一个组合valfill向量的单个向量,然后以某种方式重复该组合向量,然后从结果中构建一个新矩阵。

重复步骤可以使用一次调用来完成,rep()因为它支持向量化重复计数。换句话说,对于给定的向量输入x,它可以采用一个向量输入,times该向量输入指定重复 的每个元素的次数x。因此,挑战就变成了构建适当的x论点times

因此,解决方案首先提取 的valfillX

X[,c('val','fill')];
##      val fill
## [1,]  10  0.1
## [2,]  20  0.2
## [3,]  30  0.3
## [4,]  40  0.4
## [5,]  50  0.5
## [6,]  60  0.6

如您所见,由于我们已经索引了两列,所以我们仍然有一个矩阵,即使我们没有指定drop=F索引操作(参见R:提取或替换对象的部分)。正如将要看到的,这很方便。

在 R 中,矩阵的“矩阵角色”下面实际上只是一个普通的旧原子向量,矩阵的“向量角色”可以用于向量化操作。这就是我们如何将valandfill数据传递给rep()并适当地重复这些元素的方式。

但是,在执行此操作时,重要的是要准确了解如何将矩阵视为向量。答案是向量是由跨行的元素形成的,然后才的。(对于更高维的数组,随后的维数会跟随。IOW,向量的顺序是跨行,然后是列,然后是 z 切片等)

如果您仔细查看上面的矩阵,您会发现它不能用作我们的x参数rep(),因为val首先会跟随 s,然后是fills。实际上,我们可以相当容易地构造一个times参数来重复每个元素正确的次数,但结果向量将完全无序,并且无法将其重塑为所需的矩阵Y

实际上,在继续解释之前,我为什么不快速演示一下:

rep(X[,c('val','fill')],times=c(X[,'rep'],YC-X[,'rep']))
##  [1] 10.0 20.0 20.0 30.0 30.0 30.0 40.0 40.0 40.0 40.0 60.0  0.1  0.1  0.1  0.2  0.2  0.3  0.5  0.5  0.5  0.5  0.6  0.6  0.6

尽管上面的向量在所有正确的重复中都有所有正确的元素,但顺序是这样的,它不能形成所需的输出矩阵Y

所以,我们可以通过首先转置提取来解决这个问题:

t(X[,c('val','fill')]);
##      [,1] [,2] [,3] [,4] [,5] [,6]
## val  10.0 20.0 30.0 40.0 50.0 60.0
## fill  0.1  0.2  0.3  0.4  0.5  0.6

现在我们将valandfill向量相互交错,这样,当展平为向量时,当我们将它作为参数传递给内部使用它作为向量的函数时,就会发生这种情况,例如我们将使用rep()'sx参数,我们将以正确的顺序获取val相应fill的值,以便从中重建矩阵。让我通过将矩阵显式展平为一个向量来展示它的外观(如您所见,这种“展平”可以通过一个简单的c()调用来完成):

c(t(X[,c('val','fill')]));
##  [1] 10.0  0.1 20.0  0.2 30.0  0.3 40.0  0.4 50.0  0.5 60.0  0.6

所以,我们有我们的x论点。现在我们只需要构造times参数。

这实际上是相当棘手的。首先,我们可以认识到val值的重复计数直接在 的rep列中提供X,因此我们在X[,'rep']. 并且这些值的重复计数可以根据我在 中捕获的输出矩阵中的列数与前面提到的或 IOW的重复计数fill之间的差异来计算。问题是,我们需要交错这两个向量以符合我们的论点。YYCvalYC-X[,'rep']x

我不知道在 R 中交错两个向量的任何“内置”方式;似乎没有任何功能可以做到这一点。在处理这个问题时,我为这项任务提出了两种不同的可能解决方案,其中一种在性能和简洁性方面似乎都更好。但是由于我编写了我的原始解决方案来使用“更差”的解决方案,并且只是后来(实际上在写这个解释时)才想到第二种和“更好”的解决方案,我将在这里解释这两种方法,从第一种和更糟糕的开始一。

交错解决方案#1

交错两个向量可以通过顺序组合向量来完成,然后用精心设计的索引向量索引该组合向量,该索引向量基本上从组合向量的前半部分来回跳跃到后半部分,依次拉出每个元素每一半以交替的方式。

为了构造这个索引向量,我从一个长度等于组合向量长度一半的顺序向量开始,每个元素重复一次:

rep(1:nrow(X),each=2);
##  [1] 1 1 2 2 3 3 4 4 5 5 6 6

接下来,我添加一个二元素向量,该向量由0组合向量的一半长度组成:

nrow(X)*0:1;
## [1] 0 6

第二个加数循环通过第一个加数,实现我们需要的交错:

rep(1:nrow(X),each=2)+nrow(X)*0:1;
##  [1]  1  7  2  8  3  9  4 10  5 11  6 12

因此我们可以索引组合的重复向量来得到我们的times论点:

c(X[,'rep'],YC-X[,'rep'])[rep(1:nrow(X),each=2)+nrow(X)*0:1];
##  [1] 1 3 2 2 3 1 4 0 0 4 1 3

交错解决方案#2

交错两个向量也可以通过将这两个向量组合成一个矩阵然后再次展平它们来完成,以使它们自然地成为交错的。我相信最简单的方法是将rbind()它们放在一起,然后立即将它们压平c()

c(rbind(X[,'rep'],YC-X[,'rep']));
##  [1] 1 3 2 2 3 1 4 0 0 4 1 3

根据一些粗略的性能测试,似乎解决方案 #2 的性能更高,并且可以清楚地看到它更简洁。此外,可以很容易地将其他向量添加到rbind()调用中,但是添加到解决方案 #1 会涉及更多内容(几个增量)。

性能测试(使用大型数据集):

il1 <- function() c(X[,'rep'],YC-X[,'rep'])[rep(1:nrow(X),each=2)+nrow(X)*0:1];
il2 <- function() c(rbind(X[,'rep'],YC-X[,'rep']));
identical(il1(),il2());
## [1] TRUE
system.time({ replicate(30,il1()); });
##    user  system elapsed
##   3.750   0.000   3.761
system.time({ replicate(30,il1()); });
##    user  system elapsed
##   3.810   0.000   3.815
system.time({ replicate(30,il2()); });
##    user  system elapsed
##   1.516   0.000   1.512
system.time({ replicate(30,il2()); });
##    user  system elapsed
##   1.500   0.000   1.503

所以完整的rep()调用以正确的顺序为我们提供了我们的数据:

rep(t(X[,c('val','fill')]),times=c(rbind(X[,'rep'],YC-X[,'rep'])));
##  [1] 10.0  0.1  0.1  0.1 20.0 20.0  0.2  0.2 30.0 30.0 30.0  0.3 40.0 40.0 40.0 40.0  0.5  0.5  0.5  0.5 60.0  0.6  0.6  0.6

最后一步是使用 构建一个矩阵byrow=T,因为这就是数据最终从rep(). 而且我们还必须指定所需的行数,这与输入矩阵相同,XR(或者,如果需要,我们可以指定列数YC,或者甚至两者都指定):

Y <- matrix(rep(t(X[,c('val','fill')]),times=c(rbind(X[,'rep'],YC-X[,'rep']))),XR,byrow=T);
Y;
##      [,1] [,2] [,3] [,4]
## [1,] 10.0  0.1  0.1  0.1
## [2,] 20.0 20.0  0.2  0.2
## [3,] 30.0 30.0 30.0  0.3
## [4,] 40.0 40.0 40.0 40.0
## [5,]  0.5  0.5  0.5  0.5
## [6,] 60.0  0.6  0.6  0.6

我们完成了!

于 2015-05-12T01:35:01.657 回答