Mathematical tools for natural sciences

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.

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

The rows and columns of a matrix can be named using functions

> 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

To get the dimensions of a matrix, use the

> dim(amat)

[1] 3 4

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

> nrow(amat)

[1] 3

> ncol(amat)

[1] 4

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

[1] 5.5 6.4 8.6 9.2

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

[1] 5.5 6.4 8.6 9.2

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

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

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

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

The function call

> 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

Given a vector 'X', the function call

> 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

For a square matrix 'A', the function call

> 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

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

> 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

Given a matrix 'A' and a vector 'b' that follow the matrix equation $AX = b$, the function call

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

[1] 1 2 3

For a square matrix 'A', the function call

> 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

Similarly, the summation of every row or column can be obtained as a vector using

> 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)

[1] 5.033333 4.566667 6.933333 4.066667

> > colMeans(A)

[1] 5.100 3.875 6.475

> > rowSums(A)

[1] 15.1 13.7 20.8 12.2

> > colSums(A)

[1] 20.4 15.5 25.9

The eigenvalues and eigenvectors of a square matrix are returned by the function

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

[1] 33.148614 -4.906912 -4.000000 1.758298

> > 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.