CountBio

Matrices

A Matrix is a 2 dimensional rectangular array of numbers, symbols or expressions. Like an array, all elements of a matrix are of same data type and all the rows should be of same length.

Creating a matrix from a vector

The general format for creating a matrix of 'r' rows and 'c' columns from a vector 'vec' is,

amat <- matrix(vec, nrow=r, ncol=c, byrow=FALSE)

Here, byrow=TRUE indicates that the matrix should be filled by rows. Similarly, byrow=FALSE will fill the matrix by columns (the default).

See the examples below:

 >  x <- c(10,20,30,40,50,60,70,80,90,100,110,120)
 >  amat <- matrix(x, nrow=3, ncol=4)
 >  amat

     [,1] [,2] [,3] [,4]
[1,]   10   40   70  100
[2,]   20   50   80  110
[3,]   30   60   90  120

> amat <- matrix(x, nrow=3, ncol=4, byrow=TRUE)
> amat

     [,1] [,2] [,3] [,4]
[1,]   10   20   30   40
[2,]   50   60   70   80
[3,]   90  100  110  120

Naming the rows and columns of a matrix

The rows and columns of a matrix can be named using functions rownames() and colnames()

 >  rownames(amat) <- c("r1","r2","r3")
 >  colnames(amat) <- c("c1", "c2", "c3", "c4")
 >  amat

   c1  c2  c3  c4
r1 10  20  30  40
r2 50  60  70  80
r3 90 100 110 120

getting the number of rows and columns in a matrix

To get the dimensions of a matrix, use the dim() function:

> dim(amat)

The number of rows and columns in matrix 'amat' will be returned by function calls nrow(amat) and ncol(amat)

 >  nrow(amat)

 >  ncol(amat)

Matrix functions

R provides many functions for matrix operations. Most of the matrix computations in linear algebra can be effortlessly performed by these functions.

In order to demonstrate the matrix functions, we will create three arbitrary matrices 'A','B' and 'C' along with two vectors 'X' and 'Y'. These matrices will be used in the demonstrations ahead:

 >  ax = c(2,5,3,7,1,8,9,10,1,12,5,10,4,17,15,11)
 >  A = matrix(ax, nrow=4, ncol=4)
 >  A

     [,1] [,2] [,3] [,4]
[1,]    2    1    1    4
[2,]    5    8   12   17
[3,]    3    9    5   15
[4,]    7   10   10   11

 >  bx = c(12,5,3,17,1,18,9,10,1,12,5,10,4,15,15,14)
 >  B = matrix(bx, nrow=4, ncol=4)
 >  B

     [,1] [,2] [,3] [,4]
[1,]   12    1    1    4
[2,]    5   18   12   15
[3,]    3    9    5   15
[4,]   17   10   10   14

 >  cx = c(13,6,10,8,4,7,31,8)
 >  C = matrix(cx, nrow=4, ncol=2)
 >  C

     [,1] [,2]
[1,]   13    4
[2,]    6    7
[3,]   10   31
[4,]    8    8

 >  X = c(5, 6, 8, 9)
 >   X

 >  Y = c(8, 10, 12, 5)
 >   X

Element-wise multiplication of 2 matrices

The two matrices A and B defined above can be multiplied element-wise using the multiplication operator '*':

 >  R = A*B
 >  R

     [,1] [,2] [,3] [,4]
[1,]   24    1    1   16
[2,]   25  144  144  255
[3,]    9   81   25  225
[4,]  119  100  100  154

Matrix multiplication

The matrix multiplication of A and B is represented by the operator "%*%" as shown below:

 >  M = A %*% B
 >  
 >  M

     [,1] [,2] [,3] [,4]
[1,]  100   69   59   94
[2,]  425  427  331  558
[3,]  351  360  286  432
[4,]  351  387  287  482

Outer product of two vectors (tensor product)

The outer product of two vectors X and Y can be performed using the oprator "%o%" :

 >  X
[1] 5 6 8 9
 >  Y
[1]  8 10 12  5
 >  
 >  P = X %o% Y
 >  
 >  P

     [,1] [,2] [,3] [,4]
[1,]   40   50   60   25
[2,]   48   60   72   30
[3,]   64   80   96   40
[4,]   72   90  108   45

The above concept can be extended to perform the inner product of two matrices A and B as A %o% B

Inner product of two vectors (dot product)

The inner product(dot product) of two vector is the component-wise multiplication and the sum of it, resulting in a number. We can perform this as a simple multiplication of two vectors, already seen before:

 >  X
[1] 5 6 8 9
 >  
 >  Y
[1]  8 10 12  5
 >  
 >  D = X*Y
 >  
 >  D

[1] 40 60 96 45

Transpose of a matrix

The function call t(A) returns a matrix which is the transpose of the matrix 'A':

 >  A
     [,1] [,2] [,3] [,4]
[1,]    2    1    1    4
[2,]    5    8   12   17
[3,]    3    9    5   15
[4,]    7   10   10   11
 >  
 >  T = t(A)
 >  
 >  T

     [,1] [,2] [,3] [,4]
[1,]    2    5    3    7
[2,]    1    8    9   10
[3,]    1   12    5   10
[4,]    4   17   15   11

Creating a Diagonal matrix from a vector

Given a vector 'X', the function call diag(X) returns a square matrix with elements of 'X' as its principal diagonal :

 >  X = c(10,20,30,40,50)
 >  
 >  D = diag(X)
 >  
 >  D

     [,1] [,2] [,3] [,4] [,5]
[1,]   10    0    0    0    0
[2,]    0   20    0    0    0
[3,]    0    0   30    0    0
[4,]    0    0    0   40    0
[5,]    0    0    0    0   50

Getting the diagonal elements of a square matrix

For a square matrix 'A', the function call diag(A) returns the diagonal elements of 'A' as a vector:

 >  A
     [,1] [,2] [,3] [,4]
[1,]    2    1    1    4
[2,]    5    8   12   17
[3,]    3    9    5   15
[4,]    7   10   10   11
 >  
 >  d = diag(A)
 >  
 >  d

[1]  2  8  5 11

Creating an identity matrix

If 'k' is a scalar, we can create an identity matrix of dimension 'nxn' by the call diag(k)

 >  k = 5
 >  
 >  ds = diag(k)
 >  
 >  ds

     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    1

Solving matrix equation

Given a matrix 'A' and a vector 'b' that follow the matrix equation $AX = b$, the function call solve(A, b) returns the solution vector 'X'.

We can use this function to solve non-homogeneous linear algebraic equations. For example, given the set of three equations,

                       3x + 4y - 2z = 5 

                       4x - 5y + z  = -3 

                       10x - 6y + 5z = 13 

we can solve them as follows:

 >  xvec = c(3, 4, -2, 4, -5, 1, 10, -6, 5)
 >  
 >  A = matrix(xvec, nrow=3, ncol=3, byrow=TRUE)
 >  
 >  A
     [,1] [,2] [,3]
[1,]    3    4   -2
[2,]    4   -5    1
[3,]   10   -6    5
 >  
 >  b = c(5, -3, 13)
 >  
 >  X = solve(A, b)
 >  
 >  X

Find the inverse of a square matrix

For a square matrix 'A', the function call solve(A) returns the inverse of 'A':

 >  A
     [,1] [,2] [,3]
[1,]    3    4   -2
[2,]    4   -5    1
[3,]   10   -6    5
 >  
 >  invA = solve(A)
 > 
 >  invA 

            [,1]        [,2]       [,3]
[1,]  0.12751678  0.05369128 0.04026846
[2,]  0.06711409 -0.23489933 0.07382550
[3,] -0.17449664 -0.38926174 0.20805369

In order to check the result, we multiply the matrix invA with A to get a unit matrix within numerical approximation. See below:

 >  invA %*% A

             [,1]          [,2]          [,3]
[1,] 1.000000e+00 -1.110223e-16 -5.551115e-17
[2,] 2.220446e-16  1.000000e+00  5.551115e-17
[3,] 4.440892e-16  0.000000e+00  1.000000e+00

The mean and sum along rows and columns

For a given matrix 'A', the mean values along every row or column can be returned as a vector using functions rowMeans() and colMeans() .

Similarly, the summation of every row or column can be obtained as a vector using rowSums() and colSums() functions:

 > 
 xx = c(3.2, 5.9, 4.5, 6.8, 2.3, 1.1, 8.2, 3.9, 9.6, 6.7, 8.1, 1.5)
 >  
 >  A = matrix(xx, nrow=4, ncol=3)
 > 
 >  A

     [,1] [,2] [,3]
[1,]  3.2  2.3  9.6
[2,]  5.9  1.1  6.7
[3,]  4.5  8.2  8.1
[4,]  6.8  3.9  1.5

 >   
 >   rowMeans(A)

 >   
 >   colMeans(A)

 >   
 >   rowSums(A)

 >   
 >   colSums(A)

Eigenvalues and eigenvectors of a square matrix

The eigenvalues and eigenvectors of a square matrix are returned by the function eigen() .

 ax = c(2,5,3,7,1,8,9,10,1,12,5,10,4,17,15,11)
 >  
 >  A = matrix(ax, nrow=4, ncol=4)
 >  
 >  A

     [,1] [,2] [,3] [,4]
[1,]    2    1    1    4
[2,]    5    8   12   17
[3,]    3    9    5   15
[4,]    7   10   10   11

 >  
 >  res = eigen(A)
 >  
 >  res$values

 >  
 >  res$vectors

          [,1]       [,2]          [,3]       [,4]
[1,] 0.1086679  0.1515296 -7.789760e-16  0.8531845
[2,] 0.6433178 -0.3308090 -7.071068e-01 -0.3732923
[3,] 0.5140232  0.8457024  7.071068e-01 -0.3414212
[4,] 0.5568784 -0.3903738  1.557952e-15  0.1271242

In the above matrix of eigen vectors, each column is an eigen vector with corresponding eigen value. Thus, the first column [,1] is the eigenvector with eigenvalue 33.148614. Similarly, the second column [,2] is an eigenvector with eigenvalue -4.90692 etc.