Efficient way to create a circulant matrix in R

Question

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

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

The second row is the same as the first row except the last element is at the beginning, and so on.

Now I have the vector, say, (1, 2, 3, 4) and I want to find a efficient (fast) way to create this matrix. In practice, the 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 generate this kind of matrices over and over. A small improvement for one step will make a big difference. Thank you.

Answer 1

Give this a shot:

u<-1:4

C<-t(embed(c(u,u), length(u)+1)[,1:length(u)])
C
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    4    1    2    3
[3,]    3    4    1    2
[4,]    2    3    4    1

I should mention that this use of embed- transpose or not- should work with any input vector u.

Answer 2

We can try embed from the stats packages

f <- function(n)  t(embed(seq.int(2 * n - 1), n) %% n + 1)

and we can see

> f(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

> f(5)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    5    1    2    3    4
[3,]    4    5    1    2    3
[4,]    3    4    5    1    2
[5,]    2    3    4    5    1

> f(6)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    6    1    2    3    4    5
[3,]    5    6    1    2    3    4
[4,]    4    5    6    1    2    3
[5,]    3    4    5    6    1    2
[6,]    2    3    4    5    6    1

Answer 3

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 that Ndoogan.Circulant(1:100) takes on my machine.

Answer 4

Here are some benchmarks of 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 benchmark

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))    

Answer 5

This makes use of vector recycling (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" 

Answer 6

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

Answer 7

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

Might be faster inside a function if you constructed the double-length vector outside the sapply loop.