R / Index lower triangle vector by pairwise indices

2
votes

I am looking for a way to index a lower triangle vector of a symmetric matrix by pairwise indices in the language R. Here is some explanatory code:

Background

M is an n times n pairwise matrix, here with random data (just to show):

n <- 5
set.seed(0815)
M <- matrix(sample(n^2),n,n)
M[upper.tri(M)] <- t(M)[upper.tri(M)]
diag(M) <- 0
M

#      [,1] [,2] [,3] [,4] [,5]
# [1,]    0    6    3   18   10
# [2,]    6    0   11   23    9
# [3,]    3   11    0    5   21
# [4,]   18   23    5    0   16
# [5,]   10    9   21   16    0

Data

My data (x) is the lower triangle vector from the matrix M:

x <- M[lower.tri(M)]

And the index pairs I want to extract is stored in a vector named i:

i <- c(2:3,5)

Aim

My aim is now to extract the pairs from the lower triangle vector as in the following matrix example. 

aim <- M[i,i][lower.tri(M[i,i])]
aim
# [1] 11  9 21

My Solution

Since in my case M is not available and I don't want to generate it from x (memory issue), I created the following indexing function f:

f <- function (n,i) {
  r <- combn(i,2)[1,]
  c <- combn(i,2)[2,]
  n * r + c - (r * (r + 1))/2 - n
}

The Result

res <- x[f(n,i)]
identical(res, aim)
# [1] TRUE

Finally, my question:

Is there a more elegant or already builtin version in R of the function f, maybe one that also does not need n as argument, calculating it from the length of x?

Thank you in advance.

Sven

2
rmatrixindexinggeometrysymmetric
You might want to look for solution to extracting particular entries within r distance "matrices" since they are stored as the values in lower triangular order. (I'm not aware of a C-implemented function, but you might think about restricting your search to answers that also included Rcpp strategies.) – IRTFM

2 Answers

2
votes

I feel like your best best is to use something from the Matrix library to store your original matrix. It has several nice features which I expect could make dealing with matrices like this quite easier.

Here I use a dtpMatrix (dense, triangular, packed, matrix), without the diagonal. This takes the same amount of space in memory as your x plus just a little more for keeping track of the matrix size, etc. By using this you can refer to the rows and columns in the usual way. I then can take the desired rows and columns from your i, note though that the -1 is needed because I'm not storing the diagonal. Also notice that the @x returns just the lower triangular part of the matrix, because that's how it's stored internally.

library(Matrix)
k <- length(x)
n <- as.integer(((sqrt(8 * k + 1) + 1) / 2) - 1)
M2 <- new("dtpMatrix", uplo="L", diag="N", Dim=c(n, n), x=x)
M2[i[-1] - 1, i[-length(i)]]@x
## [1] 11  9 21

As an aside, note that this also calculates the n from the length of x, as mentioned.

0
votes

After looking at this again I'm puzzled that you didn't just do:

x[i]
#[1]  3 18 11

If you didn't want the intermediate value assigned to x you could just use:

M[lower.tri(M)][i]
#[1]  3 18 11