Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

The diag() function in R

Tags:

r

Is there a way to use the diag() function in a Matrix without using the built-in function or iteration?

   M<-matrix(1:9, ncol=3) # make a matrix 

    q5b<-function(M){ #function

    }

I know that M[1,1], M[2,2], and M[3,3] will give me the same output as diag(M). However, I can't think of a way to do this without a for loop.

My thought process was I should have a condition where row index == column index in the Matrix then print that value. I appreciate any suggestions.

like image 676
user3018479 Avatar asked Dec 10 '13 08:12

user3018479


2 Answers

You can use the functions row and col to find the indices where the column number is identical to the row number:

row(M) == col(M)
#       [,1]  [,2]  [,3]
# [1,]  TRUE FALSE FALSE
# [2,] FALSE  TRUE FALSE
# [3,] FALSE FALSE  TRUE

M[row(M) == col(M)]
# [1] 1 5 9
like image 60
Sven Hohenstein Avatar answered Sep 29 '22 00:09

Sven Hohenstein


Just subset based on another matrix:

> diag(M)
[1] 1 5 9
> M[matrix(rep(sequence(ncol(M)), 2), ncol = 2)]
[1] 1 5 9

The above would run into a problem in a non-square matrix, so we modify it as below.

As your function, one answer for question 5b could be:

q5b <- function(M) { 
  A <- sequence(ncol(M))[sequence(min(nrow(M), ncol(M)))]
  M[cbind(A, A)]
}

Update: Benchmarks are always fun

library(microbenchmark)

fun1 <- function(M) diag(M)
fun2 <- function(M) M[row(M) == col(M)]
fun3 <- function(M) {
  A <- sequence(ncol(M))[sequence(min(nrow(M), ncol(M)))]
  M[cbind(A, A)]
}    

set.seed(1)
M <- matrix(rnorm(1000*1000), ncol = 1000)

microbenchmark(fun1(M), fun2(M), fun3(M), times = 100)
# Unit: microseconds
#     expr       min        lq     median        uq        max neval
#  fun1(M)  4654.825  4747.408  4822.8865  4912.690   5877.866   100
#  fun2(M) 53270.266 54813.606 55059.0695 55749.062 200384.531   100
#  fun3(M)    66.284    82.321   118.8835   129.361    191.155   100
like image 36
A5C1D2H2I1M1N2O1R2T1 Avatar answered Sep 28 '22 23:09

A5C1D2H2I1M1N2O1R2T1