如何在向量中创建具有不同重复值的矩阵

How to create a matrix with different repeats of values in a vector

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

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

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

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

我的目标是生成一个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,但还是不行。

输出矩阵的每一行都可以通过调用 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 矩阵,希望它适合您的系统内存。但是,如果 运行 系统内存不足,您可能无法将它推得太大。

解决方案

这并不容易,但我想出了一种方法来完成此任务,使用对 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 fill 分配给 X,而不是问题中指定的 a b c,并且我在解决方案中使用了这些列名索引 X(而不是使用数字索引),因为我通常更喜欢尽可能提高人类可读性,但就解决方案的正确性和性能而言,这个细节可以忽略不计。)

性能

这实际上比@josilber 的解决方案具有显着的性能优势,因为他使用 apply() 在内部循环矩阵的行(在 R 语言中传统上称为 "hidden loop"),而我的解决方案的核心是对 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 的每个元素多少次。因此,挑战就变成了构建适当的 xtimes 参数。

因此,解决方案首先提取 Xvalfill 列:

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: Extract or Replace Parts of an Object)。这很方便,稍后会看到。

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

但是,在执行此操作时,重要的是要准确理解如何将矩阵视为向量。答案是向量由以下元素组成 跨行 并且只有 跨列 。 (对于高维数组,随后是后续维度。IOW,向量的顺序是跨行,然后是列,然后是 z 切片等)

如果你仔细看上面的矩阵,你会发现它不能用作我们 rep()x 参数,因为 val 将首先出现,然后是 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

现在我们有 valfill 向量相互交错,这样,当扁平化为向量时,当我们将它作为参数传递给内部函数时会发生这种情况将它用作向量,例如我们将使用 rep()x 参数,我们将以正确的顺序获得 val 和相应的 fill 值以进行重建它们中的一个矩阵。让我通过将矩阵显式展平为向量来演示这一点(如您所见,这个 "flattening" 可以通过简单的 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 值的重复计数直接在 Xrep 列中提供,所以我们在 X[,'rep'] 中有它。 fill 值的重复计数可以根据我在 YC 中捕获的输出矩阵 Y 中的列数与上述重复计数之间的差异来计算对于 val,或 IOW,YC-X[,'rep']。问题是,我们需要交错这两个向量以符合我们的 x 论点。

我不知道 "built-in" 在 R 中交错两个向量的方法;似乎没有任何功能可以做到这一点。在解决这个问题时,我为这个任务提出了两种不同的可能解决方案,其中一种似乎在性能和简洁性方面都更好。但是因为我写了我原来的解决方案来使用 "worse" 一个,只是后来(实际上是在写这个解释时)想到了第二个和 "better" 一个,我将在这里解释这两种方法,从第一个也是更糟糕的一个。

交织解决方案 #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

大功告成!