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.

Selsyn

If you put your code into a function, you can compile it to bytecode with the compiler library (shipped with R), and this should give you some speed gain. for more details: dirk.eddelbuettel.com/blog/2011/04/12 – Perianth

And if this is going to be used repeatedly, using package Rcpp would greatly improve performance. – Baleful

Answer 1

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

Utah

Answer 2

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" 

Pete

Many times slower than the OP's solution on my machine. – Baleful

@DWin As others are seeing differences, I observe that mine is only second to yours. Thanks again. – Pete

I thought I was out of the running. The OP's solution beats mine on my machine (with my version) – Baleful

Very nice. This is what I originally tried for, but couldn't figure it out. Adding a sacrificial column to the matrix is the key insight. Thanks! – Cheddar

Answer 3

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

Ramage

@DWin -- Thanks. It could be made faster by avoiding constructing the matrix byrow, but would be at the expense of a little clarity. (Not that this isn't plenty obfuscated already!) – Cheddar

Clarity occurred when I looked at these separately: (1:4 - rep(1:4, each=4)) and rep(1:4, each=4) – Baleful

It can be made (much?) faster by working with integers, not numerics. Try replacing n + 1 with n + 1L. – Rompers

@Rompers -- Thanks for suggesting that one-letter edit. Nice catch! (Makes it about 15% faster in my timings.) – Cheddar

Answer 4

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.

Baleful

Answer 5

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.

Cleanup

nice! is the advantage similar for smaller matrices? – Pete

Answer 6

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.

Factory

Your result does not match the OP's desired matrix, please consider revising your solution so that it matches the OP's requested output. Thanks. – Penna

Ah, fair enough. A transpose should take care of that problem. Thanks! – Factory

Answer 7

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

Violative