Geek Answers Handbook

Efficient way to create a circulant matrix in R

I want to create a circulant matrix from a vector to R. A circulant matrix is ​​a matrix of the following form.

1 2 3 4 4 1 2 3 3 4 1 2 2 3 4 1

The second line is the same as the first line, except that the last element is at the beginning, etc.

Now I have a vector, say (1, 2, 3, 4), and I want to find an efficient (fast) way to create this matrix. In practice, numbers are not integers and can be any numbers.

Here is what I am doing now.

 x <- 1:4 n <- length(x) mat <- matrix(NA, n, n) for (i in 1:n) { mat[i, ] <- c(x[-(1:(n+1-i))], x[1:(n+1-i)]) }

I wonder if there is a faster way to do this? I need to create these types of matrices over and over again. A slight improvement by one step will go a long way. Thanks.

+6
source share
5 answers

Here are some guidelines for suggested solutions.

ndoogan takes the lead!

Benchmark

 x <- 1:100 microbenchmark( OP.Circulant(x), Josh.Circulant(x), Dwin.Circulant(x) , Matt.Circulant(x), Matt.Circulant2(x), Ndoogan.Circulant(x), times=100 ) # Unit: microseconds # expr min lq median uq max # 1 Dwin.Circulant(x) 1232.775 1288.1590 1358.999 1504.4490 2900.430 # 2 Josh.Circulant(x) 1081.080 1086.3470 1097.863 1125.8745 2526.237 # 3 Matt.Circulant(x) 61924.920 64579.3735 65948.152 129359.7895 137371.570 # 4 Matt.Circulant2(x) 12746.096 13499.0580 13832.939 14346.8570 16308.040 # 5 Ndoogan.Circulant(x) 469.502 487.2285 528.591 585.8275 1522.363 # 6 OP.Circulant(x) 1291.352 1363.8395 1421.509 1513.4950 2714.707

Code used for testing

 OP.Circulant <- function(x) { n <- length(x) mat <- matrix(NA, n, n) for (i in 1:n) { mat[i, ] <- c(x[-(1:(n + 1 - i))], x[1:(n + 1 - i)]) } return(mat) } rotn <- function(x, n) rep(x, 2)[n:(n + length(x) - 1)] Dwin.Circulant <- function(x) { n <- length(x) return(t(sapply(x[c(1L, n:2)], rotn, x = x))) } Josh.Circulant <- function(x, nrow = length(x)) { m <- length(x) return(matrix(x[(1:m - rep(1:nrow, each = m))%%m + 1L], ncol = m, byrow = TRUE)) } Matt.Circulant <- function(x) { n <- length(x) mat <- matrix(, n, n) for (i in seq(-n + 1, n - 1)) { mat[row(mat) == col(mat) - i] = x[i%%n + 1] } return(mat) } Matt.Circulant2 <- function(x) { n <- length(x) return(rbind(x[], do.call(rbind, lapply(seq(n - 1), function(i) c(tail(x, i), head(x, -i)))))) } Ndoogan.Circulant <-function(x) { n <- length(x) suppressWarnings( matrix(x[matrix(1:n,n+1,n+1,byrow=T)[c(1,n:2),1:n]],n,n)) } # check for identical results (all TRUE) check <- OP.Circulant(x) identical(check, OP.Circulant(x)) identical(check, Dwin.Circulant(x)) identical(check, Josh.Circulant(x)) identical(check, Matt.Circulant(x)) identical(check, Matt.Circulant2(x)) identical(check, Ndoogan.Circulant(x))
+4
source

This uses vector processing (it throws a warning):

 circ<-function(x) { n<-length(x) matrix(x[matrix(1:n,n+1,n+1,byrow=T)[c(1,n:2),1:n]],n,n) } circ(letters[1:4]) # [,1] [,2] [,3] [,4] #[1,] "a" "b" "c" "d" #[2,] "d" "a" "b" "c" #[3,] "c" "d" "a" "b" #[4,] "b" "c" "d" "a"
+5
source 
 rotn <- function(x,n) rep(x,2)[n:(n+length(x)-1)] sapply(c(1,4:2), rotn, x=1:4) [,1] [,2] [,3] [,4] [1,] 1 4 3 2 [2,] 2 1 4 3 [3,] 3 2 1 4 [4,] 4 3 2 1

Be faster inside the function if you built a double-length vector outside the sapply loop.

+4
source 
 circulant <- function(x, nrow = length(x)) { n <- length(x) matrix(x[(1:n - rep(1:nrow, each=n)) %% n + 1L], ncol=n, byrow=TRUE) } circulant(1:4) # [,1] [,2] [,3] [,4] # [1,] 1 2 3 4 # [2,] 4 1 2 3 # [3,] 3 4 1 2 # [4,] 2 3 4 1 circulant(7:9, nrow=5) # [,1] [,2] [,3] # [1,] 7 8 9 # [2,] 9 7 8 # [3,] 8 9 7 # [4,] 7 8 9 # [5,] 9 7 8 circulant(10:1, nrow=2) # [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] # [1,] 10 9 8 7 6 5 4 3 2 1 # [2,] 1 10 9 8 7 6 5 4 3 2
+4
source

Here is a solution using Rcpp:

 library(Rcpp) cppFunction(" IntegerMatrix myCirculant(const int n) { IntegerMatrix res(n); int val = 1; int dval = 2; for (int i = 0; i < n*n; i++) { res[i] = val; if (val > 1) { if (val != dval) { val--; } else { if (dval == n) { dval = 1; } else { dval++; } } } else { val = n; } } return res; }") myCirculant(100)

works only for integers and takes 1/10 of the time when Ndoogan.Circulant(1:100) takes my car.

0
source

More articles:


All Articles