13

I've got a large nxn matrix and would like to take off-diagonal slices of varying sizes. For example:

1 2 3 4 5 6
1 2 3 4 5 6
1 2 3 4 5 6
1 2 3 4 5 6
1 2 3 4 5 6
1 2 3 4 5 6

I'd like an R function which, when given the matrix and "width of diagonal slice" would return an nxn matrix of just those values. So for the matrix above and, say, 3, I'd get:

1 x x x x x
1 2 x x x x
1 2 3 x x x
x 2 3 4 x x
x x 3 4 5 x
x x x 4 5 6

At the moment I'm using (forgive me) a for loop which is incredibly slow:

getDiags<-function(ndiags, cormat){
  resmat=matrix(ncol=ncol(cormat),nrow=nrow(cormat))
  dimnames(resmat)<-dimnames(cormat)
  for(j in 1:ndiags){
    resmat[row(resmat) == col(resmat) + j] <- 
      cormat[row(cormat) == col(cormat) + j]
  }
  return(resmat)
}

I realise that this is a very "un-R" way to go about solving this problem. Is there a better way to do it, probably using diag or lower.tri?

blmoore
  • 1,487
  • 16
  • 31

3 Answers3

15
size <- 6
mat <- matrix(seq_len(size ^ 2), ncol = size)


low <- 0
high <- 3

delta <- rep(seq_len(ncol(mat)), nrow(mat)) - 
    rep(seq_len(nrow(mat)), each = ncol(mat))
#or Ben Bolker's better alternative
delta <- row(mat) - col(mat)
mat[delta < low | delta > high] <- NA
mat

this works with 5000 x 5000 matrices on my machine

Thierry
  • 18,049
  • 5
  • 48
  • 66
  • This method for 10 diags in my matrix: 2.125 My old method (elapsed): 15.170 - Great, thanks! – blmoore Aug 01 '12 at 13:33
  • 10
    You can also use `delta <- row(mat)-col(mat)`; it's probably a little slower to run, but is very easy to read ... – Ben Bolker Aug 01 '12 at 15:48
  • 1
    @BenBolker +1 for your solution. It is actually about 30% faster and it can handle larger matrices (8000 vs 6000 on my machine) – Thierry Aug 02 '12 at 07:53
1

If you want to use upper.tri and lower.tri you could write functions like these:

cormat <- mapply(rep, 1:6, 6)

u.diags <- function(X, n) {
  X[n:nrow(X),][lower.tri(X[n:nrow(X),])] <- NA
  return(X)
}

or

l.diags <- function(X, n) {
  X[,n:ncol(X)][upper.tri(X[,n:ncol(X)])] <- NA
  return(X)
}

or

n.diags <- function(X, n.u, n.l) {
  X[n.u:nrow(X),][lower.tri(X[n.u:nrow(X),])] <- NA
  X[,n.l:ncol(X)][upper.tri(X[,n.l:ncol(X)])] <- NA
  return(X)
}

l.diags(cormat, 3)
u.diags(cormat, 3)
n.diags(cormat, 3, 1)
Bryce F
  • 101
  • 1
  • 1
  • 5
0

you can do:

matrix:

m<-
matrix(1:6,ncol = 6, nrow=6 ,byrow = T)

function:

n_diag <- function (x, n) {
    d <- dim(x)
    ndiag <- .row(d) - n >= .col(d)
    x[upper.tri(x) | ndiag] <- NA
    return(x)
}

call:

n_diag(m,3)

#     [,1] [,2] [,3] [,4] [,5] [,6]
#[1,]    1   NA   NA   NA   NA   NA
#[2,]    1    2   NA   NA   NA   NA
#[3,]    1    2    3   NA   NA   NA
#[4,]   NA    2    3    4   NA   NA
#[5,]   NA   NA    3    4    5   NA
#[6,]   NA   NA   NA    4    5    6

just for fun:

#lapply(1:6, n_diag, x = m)
Andre Elrico
  • 10,956
  • 6
  • 50
  • 69