【问题标题】:What is the fast way to turn this vector into symmetric matrix?将此向量转换为对称矩阵的快速方法是什么?
【发布时间】:2020-11-05 02:29:10
【问题描述】:

大家!假设我有一个长度为 n(n+1)/2 的向量:

a = (a_11, a_12, a_22, ...., a_nn) 

现在我想把它变成一个对称矩阵,意思是

我可以一个一个地赋值,但我想知道是否有更快的方法来创建这个矩阵?非常感谢!!

【问题讨论】:

  • 提供一个实际的向量和相应的预期输出而不是索引会让这个问题更容易回答。
  • 输入向量中的a_12是不是应该同时插入到矩阵中的a_12和a_21中?

标签: r


【解决方案1】:

你可以写一个这样的函数:

如果你的向量是a11,a12,a13..a1n,a22,a23..a2n, a33,..a3n,..ann 你可以这样做:

vec2mat <- function(x){
  p <- sqrt(1 + 8 * length(x))/ 2 - 0.5
  m <- matrix(0, p, p)
  m[lower.tri(m, diag = TRUE)] <- x
  m[upper.tri(m)] <- (t(m))[upper.tri(m)]
  m
}

现在:

vec2mat(1:6)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    4    5
[3,]    3    5    6
vec2mat(1:10)
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    2    5    6    7
[3,]    3    6    8    9
[4,]    4    7    9   10

如果你有a11, a12,a22,a31, a32, a33...

vec2mat <- function(x){
  p <- sqrt(1 + 8 * length(x))/ 2 - 0.5
  m <- matrix(0, p, p)
  m[upper.tri(m, diag = TRUE)] <- x
  m[lower.tri(m)] <- t(m)[lower.tri(m)]
  m
}
 vec2mat(1:10)
     [,1] [,2] [,3] [,4]
[1,]    1    2    4    7
[2,]    2    3    5    8
[3,]    4    5    6    9
[4,]    7    8    9   10

【讨论】:

  • 对于 1000 个元素它崩溃了......这是为什么?
  • @ManosPapadakis 因为你不能有 1000 个元素。请注意,数字必须是三角形数字。即包括对角线的方阵的下三角形或上三角形。因此,如果它是 44 行和 44 列的正方形,数字即向量的长度应为 990,如果它是 45 行和 45 列的正方形,则长度应为 1035。
  • 非常感谢!
【解决方案2】:

你可以试试:

#define n
n <- 2

#Create n(n+1)/2 objects in global environment
#OP already has that
a_11 <- 5
a_12 <- 9
a_22 <- 8

#Create n X n matrix with NA
mat <- matrix(nrow = n, ncol = n)
#Get all the individual objects in one vector
vec <- unlist(mget(ls(pattern = 'a_')), use.names = FALSE)
#Replace upper (or lower) triangular elements with it
mat[upper.tri(mat, diag = TRUE)] <- vec
#Copy the elements to other half.
mat[lower.tri(mat)] <- mat[upper.tri(mat)]

#      [,1] [,2]
#[1,]    5    9
#[2,]    9    8

【讨论】:

    【解决方案3】:

    你可以像这样创建一个稀疏矩阵:

    a <- 1:6
    n <- as.integer(-0.5 + sqrt(0.25 + 2 * length(a)))
    
    library(Matrix)
    sparseMatrix(x = a, dims = c(n, n), symmetric = TRUE, 
                 i = sequence(1:n), j = rep(1:n, 1:n))
    #3 x 3 sparse Matrix of class "dsCMatrix"
    #          
    #[1,] 1 2 4
    #[2,] 2 3 5
    #[3,] 4 5 6
    

    如果您需要密集矩阵,请在结果上使用as.matrix。如果您的向量的顺序与您显示的不同(例如,其他一些答案假设的主要行),您需要稍微调整 ij 的计算。

    【讨论】:

      【解决方案4】:
      a <- 1:25
      sq.m <- matrix( a, ncol=sqrt(length(a) ) )
      
      > sq.m
           [,1] [,2] [,3] [,4] [,5]
      [1,]    1    6   11   16   21
      [2,]    2    7   12   17   22
      [3,]    3    8   13   18   23
      [4,]    4    9   14   19   24
      [5,]    5   10   15   20   25
      

      显然,这不是一个对称矩阵,但如果初始向量的顺序是一个,那么它就会成功。

      如果您想强制一个非方阵,使其上三角元素与下三角元素相同,则可以在赋值的两侧进行索引:

       sq.m[ upper.tri(sq.m) ] <- sq.m[lower.tri(sq.m)]
      > sq.m
           [,1] [,2] [,3] [,4] [,5]
      [1,]    1    2    3    5   10
      [2,]    2    7    4    8   14
      [3,]    3    8   13    9   15
      [4,]    4    9   14   19   20
      [5,]    5   10   15   20   25
      

      【讨论】:

      • OP 在输入向量中没有 n^2 个元素。我认为目的是在对角线上复制某些元素以创建一个
      • 他的问题演示告诉我他确实有 n^2 个元素。你有什么相反的证据?
      • Suppose I have a vector of length n(n+1)/2:
      • 所以你只有下三角元素?啊,现在我明白了。他只有对角线和上三角元素。
      • 没关系,是对称矩阵。可以是上三角也可以是下三角。
      【解决方案5】:

      如果你有大向量并且速度很重要,那么你可以使用这个使用 Rfast 包的函数(取自@Onyambu)。

      vec2mat_rfast <- function(x){
          p <- sqrt(1 + 8 * length(x))/ 2 - 0.5
          m <- matrix(0, p, p)
          m[Rfast::upper_tri(m, diag = TRUE)] <- x
          m[Rfast::lower_tri(m)] <- Rfast::transpose(m)[Rfast::lower_tri(m)]
          m
      }
      
      vec2mat <- function(x){
          p <- sqrt(1 + 8 * length(x))/ 2 - 0.5
          m <- matrix(0, p, p)
          m[upper.tri(m, diag = TRUE)] <- x
          m[lower.tri(m)] <- t(m)[lower.tri(m)]
          m
      }
      
      y=numeric(500500)
      
      
      > microbenchmark::microbenchmark(a<-vec2mat(y),b<-vec2mat_rfast(y),times = 10)
      Unit: milliseconds
                          expr     min       lq      mean    median       uq      max neval
               a <- vec2mat(y) 88.9461 101.5294 114.96752 108.86680 112.9854 180.9679    10
         b <- vec2mat_rfast(y) 33.1351  34.5294  49.61295  45.26955  62.8214  80.5069    10
      
      > all.equal(a,b)
      [1] TRUE
      

      【讨论】:

        猜你喜欢
        • 1970-01-01
        • 2019-01-15
        • 2021-01-02
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 2020-09-02
        • 2020-10-01
        相关资源
        最近更新 更多