CODEWITHSUNDEEP
r
0o0o0o0
0o0o0o0

Reputation: 35

What is the fast way to turn this vector into symmetric matrix?

everyone! Suppose I have a vector of length n(n+1)/2:

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

Now I'd like to turn it into a symmetric matrix, mean

enter image description here

I could assign the value one by one, but I'm wondering if there is some faster to create this matrix? Thanks so much!!

Upvotes: 1

Views: 903

Answers (6)

Frank Harrell
Frank Harrell

Reputation: 2230

It would be interesting to compare the speed to the approaches above with this:

# For i,j element in pxp symmetric matrix y compute subscript in compact vector x (columns moving fastest)
ssub <- function(i, j, p) {
  ii <- pmin(i, j)
  jj <- pmax(i, j)
  (ii - 1) * p + jj - (ii - 1) - (ii - 1) * (ii - 2) / 2
}

p <- 5
y <- matrix(NA, p, p)
nv <- p * (p + 1) / 2
x <- 1 : nv
y[] <- x[ssub(row(y), col(y), p)]
y

Upvotes: 0

Manos Papadakis
Manos Papadakis

Reputation: 593

if you have big vectors and speed matters then you can use this function (taken from @Onyambu) that uses Rfast package.

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

Upvotes: 0

Roland
Roland

Reputation: 132706

You could create a sparse matrix like this:

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

Use as.matrix on the result if you need a dense matrix. If the order of your vector is different than you show (e.g., row-major as some of the other answers assume), you need to adjust calculation of i and j slightly.

Upvotes: 1

Onyambu
Onyambu

Reputation: 79208

You could write a function that does this:

IF you had your vector as a11,a12,a13..a1n,a22,a23..a2n, a33,..a3n,..ann They you could do:

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
}

Now:

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

if you had 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

Upvotes: 3

IRTFM
IRTFM

Reputation: 263342

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

Obviously this is not a symmetric matrix but if the order of the initial vector would have made one, then it would have succeeded.

If you wanted to coerce a non-square matrix so its upper triangular elements were the same as the lower triangular elements then indexing can be done on both sides of an assignment:

 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

Upvotes: 0

Ronak Shah
Ronak Shah

Reputation: 388982

You can try :

#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

Upvotes: 1

Related Questions