banbh
banbh

Reputation: 1525

Diagonal of a 3D array

How can one extract the "diagonal" from three-dimensional array in R? For a matrix (2D array) one can use the diag(...) function. In a similar way, given an N x N x M array, a natural operation is to convert it into an N x M matrix by taking the diagonal from each N x N slice and returning it as a matrix.

It's easy to do this using a loop, but that is not idiomatic R and is slow. Another possibility is to use slightly complex indexing (see my own answer to this question) but it is a bit hard to read. What other alternatives are there? Is there a standard R way to do this?

Upvotes: 2

Views: 2742

Answers (5)

Tibo
Tibo

Reputation: 68

Another approach is subseting the 3 dimensions array with a 2 dimensions matrix:

a <- array(1:100,dim = c(5,5,4))
ref <- cbind(1:5,1:5,rep(1:4,each= 5))
a[ref]

Output is a vector instead of a matrix. On my computer it is more efficient than apply() and you can also fill the diagonal values.

Upvotes: 0

IRTFM
IRTFM

Reputation: 263411

Although I'm not enamored of the term "3d.diagonal" for this result, it can be achieved with this simple function (up to identity modulo transpose):

arr <- array(1:27,c(3,3,3)  )
apply(arr, 3, function(x) x[row(x)==col(x)] )
# returns same value as diag.3d (arr)
     [,1] [,2] [,3]
[1,]    1   10   19
[2,]    5   14   23
[3,]    9   18   27

I think a "real diagonal" would be arr[ cbind(1:3,1:3,1:3) ]

Upvotes: 3

user31264
user31264

Reputation: 6737

Create an array and fill it by some values:

> a=array(0,c(10,10,5))
> for (i in 1:10) for (j in 1:10) for (k in 1:5) a[i,j,k]=100*i+10*j+k-111

Run the apply function:

> apply(a,3,diag)
      [,1] [,2] [,3] [,4] [,5]
 [1,]    0    1    2    3    4
 [2,]  110  111  112  113  114
 [3,]  220  221  222  223  224
 [4,]  330  331  332  333  334
 [5,]  440  441  442  443  444
 [6,]  550  551  552  553  554
 [7,]  660  661  662  663  664
 [8,]  770  771  772  773  774
 [9,]  880  881  882  883  884
[10,]  990  991  992  993  994

Upvotes: 4

eddi
eddi

Reputation: 49448

Various diagonals:

A = array(1:12, c(2, 2, 3))

apply(A, 1, diag)
#     [,1] [,2]
#[1,]    1    2
#[2,]    7    8
apply(A, 2, diag)
#     [,1] [,2]
#[1,]    1    3
#[2,]    6    8
apply(A, 3, diag)
#     [,1] [,2] [,3]
#[1,]    1    5    9
#[2,]    4    8   12

Upvotes: 3

banbh
banbh

Reputation: 1525

One possible approach is to use indexing, where the indices are a matrix with three columns. For example:

diag.3d <- function(A) {
  # Expect a N x N x M array

  stopifnot(length(dim(A)) == 3)
  n <- nrow(A)
  stopifnot(n == ncol(A))
  m <- dim(A)[3]

  IXS <- cbind(1:n, 1:n, rep(1:m, each = n))
  cn <- colnames(A)
  rn <- dimnames(A)[[3]]
  matrix(A[IXS], ncol = n, byrow = T, dimnames = list(rn, cn))
}

Although indices (in variable IXS) seem hard to read.

Upvotes: 1

Related Questions