【问题标题】:How to create a matrix with different repeats of values in a vector如何在向量中创建具有不同重复值的矩阵
【发布时间】:2015-07-22 13:58:13
【问题描述】:

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

假设我有一个名为 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,但还是不行。

【问题讨论】:

    标签: r matrix rep


    【解决方案1】:

    输出矩阵的每一行都可以通过对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 矩阵,希望它适合您的系统内存。但是,您可能无法在不耗尽系统内存的情况下将其推得更大。

    【讨论】:

    • @Ted apply(X, 1, f)X 的每一行上运行函数f。我们选择运行的函数是function(x) rep(x[2:3], c(x[1], 4-x[1])),它调用传递的行x,然后将该行的各个部分传递给rep函数,返回结果。您可能会发现 the following 作为 apply 函数的教程很有帮助。
    【解决方案2】:

    解决方案

    这并不容易,但我想出了一种方法来完成这项任务,只需对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 语言中传统上称为“隐藏循环”),而我的解决方案的核心是对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 中,矩阵的“矩阵角色”实际上只是一个普通的旧原子向量,并且矩阵的“向量角色”可以用于向量化操作。这就是我们可以将valfill 数据传递给rep() 并适当重复这些元素的方式。

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

    如果你仔细查看上面的矩阵,你会发现它不能用作rep()x 参数,因为vals 将首先出现,然后是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 值,以便从中重建矩阵。让我通过将矩阵显式展平为一个向量来展示这一点(如您所见,这种“展平”可以通过简单的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 的上述重复计数之间的差来计算,或者爱荷华州,YC-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
    

    我们完成了!

    【讨论】:

    • 您的回答给我留下了深刻的印象。非常感谢您所有精彩而详细的解释,这对我真的很有帮助。你和 josilber 都很棒 :) 非常感谢!
    猜你喜欢
    • 1970-01-01
    • 2021-08-08
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2022-06-14
    • 1970-01-01
    • 2021-09-15
    • 1970-01-01
    相关资源
    最近更新 更多