# Linear algebra in R

June 4, 2019
By

(This article was first published on R Programming – DataScience+, and kindly contributed to R-bloggers)

Are you interested in guest posting? Publish at DataScience+ via your editor (i.e., RStudio).

## Tags

In this article, you learn how to do linear algebra in R. In particular, I will discuss how to create a matrix in R, Element-wise operations in R, Basic Matrix Operations in R, How to Combine Matrices in R, Creating Means and Sums in R and Advanced Matrix Operations in R.

## Creating Matrix

Let us start by creating a Matrix in R:

```A <- matrix(data = 1:25, nrow = 5, ncol = 5, byrow = TRUE)
A
[,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15
[4,]   16   17   18   19   20
[5,]   21   22   23   24   25

```
```B <- matrix(data = 25:49, nrow = 5, ncol = 5, byrow = FALSE)
B
[,1] [,2] [,3] [,4] [,5]
[1,]   25   30   35   40   45
[2,]   26   31   36   41   46
[3,]   27   32   37   42   47
[4,]   28   33   38   43   48
[5,]   29   34   39   44   49

```

## Element-wise operations

Let us start move over to element-wise operations in R:

Firt we will make an Addition in R:

```A + B
[,1] [,2] [,3] [,4] [,5]
[1,]   26   32   38   44   50
[2,]   32   38   44   50   56
[3,]   38   44   50   56   62
[4,]   44   50   56   62   68
[5,]   50   56   62   68   74

```

Next we will make a Subtraction in R:

```A - B
[,1] [,2] [,3] [,4] [,5]
[1,]  -24  -28  -32  -36  -40
[2,]  -20  -24  -28  -32  -36
[3,]  -16  -20  -24  -28  -32
[4,]  -12  -16  -20  -24  -28
[5,]   -8  -12  -16  -20  -24

```

Now we will try a Multiplication in R:

```A * B
[,1] [,2] [,3] [,4] [,5]
[1,]   25   60  105  160  225
[2,]  156  217  288  369  460
[3,]  297  384  481  588  705
[4,]  448  561  684  817  960
[5,]  609  748  897 1056 1225

```

Furthermore we will make a Division in R:

```A / B
[,1]       [,2]       [,3]      [,4]      [,5]
[1,] 0.0400000 0.06666667 0.08571429 0.1000000 0.1111111
[2,] 0.2307692 0.22580645 0.22222222 0.2195122 0.2173913
[3,] 0.4074074 0.37500000 0.35135135 0.3333333 0.3191489
[4,] 0.5714286 0.51515152 0.47368421 0.4418605 0.4166667
[5,] 0.7241379 0.64705882 0.58974359 0.5454545 0.5102041

```

Let us try a Exponentiation in R:

```A ^ B
[,1]         [,2]         [,3]         [,4]         [,5]
[1,] 1.000000e+00 1.073742e+09 5.003155e+16 1.208926e+24 2.842171e+31
[2,] 1.705817e+20 1.577754e+26 3.245186e+32 1.330279e+39 1.000000e+46
[3,] 1.310999e+28 3.418219e+34 1.644008e+41 1.372074e+48 1.889249e+55
[4,] 5.192297e+33 4.025450e+40 5.015973e+47 9.691809e+54 2.814750e+62
[5,] 2.209833e+38 4.389056e+45 1.280518e+53 5.361603e+60 3.155444e+68

```

## Basic Matrix Operations

Now we will move over to basic matrix operations in R.

Let us begin with making a Transpose in R:

```t(A)
[,1] [,2] [,3] [,4] [,5]
[1,]    1    6   11   16   21
[2,]    2    7   12   17   22
[3,]    3    8   13   18   23
[4,]    4    9   14   19   24
[5,]    5   10   15   20   25

```
```t(B)
[,1] [,2] [,3] [,4] [,5]
[1,]   25   26   27   28   29
[2,]   30   31   32   33   34
[3,]   35   36   37   38   39
[4,]   40   41   42   43   44
[5,]   45   46   47   48   49

```

Now let us try a Multiplication in R:

```A %*% B
[,1] [,2] [,3] [,4] [,5]
[1,]  415  490  565  640  715
[2,] 1090 1290 1490 1690 1890
[3,] 1765 2090 2415 2740 3065
[4,] 2440 2890 3340 3790 4240
[5,] 3115 3690 4265 4840 5415

```
```B %*% A
[,1] [,2] [,3] [,4] [,5]
[1,] 2175 2350 2525 2700 2875
[2,] 2230 2410 2590 2770 2950
[3,] 2285 2470 2655 2840 3025
[4,] 2340 2530 2720 2910 3100
[5,] 2395 2590 2785 2980 3175

```

We can also mae a Outer Product in R:

```X <- matrix(data = 1:10, nrow = 10)
Y <- matrix(data = 21:30, nrow = 10)
X %o% Y
, , 1, 1

[,1]
[1,]   21
[2,]   42
[3,]   63
[4,]   84
[5,]  105
[6,]  126
[7,]  147
[8,]  168
[9,]  189
[10,]  210

, , 2, 1

[,1]
[1,]   22
[2,]   44
[3,]   66
[4,]   88
[5,]  110
[6,]  132
[7,]  154
[8,]  176
[9,]  198
[10,]  220
................

```
```X %*% t(Y)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]   21   22   23   24   25   26   27   28   29    30
[2,]   42   44   46   48   50   52   54   56   58    60
[3,]   63   66   69   72   75   78   81   84   87    90
[4,]   84   88   92   96  100  104  108  112  116   120
[5,]  105  110  115  120  125  130  135  140  145   150
[6,]  126  132  138  144  150  156  162  168  174   180
[7,]  147  154  161  168  175  182  189  196  203   210
[8,]  168  176  184  192  200  208  216  224  232   240
[9,]  189  198  207  216  225  234  243  252  261   270
[10,]  210  220  230  240  250  260  270  280  290   300

```

It is also possible to make Cross Products in R:

```crossprod(A)
[,1] [,2] [,3] [,4] [,5]
[1,]  855  910  965 1020 1075
[2,]  910  970 1030 1090 1150
[3,]  965 1030 1095 1160 1225
[4,] 1020 1090 1160 1230 1300
[5,] 1075 1150 1225 1300 1375

```
```t(A) %*% A
[,1] [,2] [,3] [,4] [,5]
[1,]  855  910  965 1020 1075
[2,]  910  970 1030 1090 1150
[3,]  965 1030 1095 1160 1225
[4,] 1020 1090 1160 1230 1300
[5,] 1075 1150 1225 1300 1375

```
```crossprod(A, B)
[,1] [,2] [,3] [,4] [,5]
[1,] 1535 1810 2085 2360 2635
[2,] 1670 1970 2270 2570 2870
[3,] 1805 2130 2455 2780 3105
[4,] 1940 2290 2640 2990 3340
[5,] 2075 2450 2825 3200 3575

```
```t(A) %*% B
[,1] [,2] [,3] [,4] [,5]
[1,] 1535 1810 2085 2360 2635
[2,] 1670 1970 2270 2570 2870
[3,] 1805 2130 2455 2780 3105
[4,] 1940 2290 2640 2990 3340
[5,] 2075 2450 2825 3200 3575

```

Now we will make a Diagonal in R:

```diag(A)
[1]  1  7 13 19 25

```
```diag(B)
[1] 25 31 37 43 49

```
```diag(c(1, 2, 3))
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    2    0
[3,]    0    0    3

```
```diag(5)
[,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

```
```diag(3)
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

```

## Combine Matrices

Let us move over to combining matrices in R.
First we will make a Horizontally combine in R:

```cbind(A, B)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    2    3    4    5   25   30   35   40    45
[2,]    6    7    8    9   10   26   31   36   41    46
[3,]   11   12   13   14   15   27   32   37   42    47
[4,]   16   17   18   19   20   28   33   38   43    48
[5,]   21   22   23   24   25   29   34   39   44    49

```

Now let us try a Vertically combine in R:

```rbind(A, B)
[,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15
[4,]   16   17   18   19   20
[5,]   21   22   23   24   25
[6,]   25   30   35   40   45
[7,]   26   31   36   41   46
[8,]   27   32   37   42   47
[9,]   28   33   38   43   48
[10,]   29   34   39   44   49

```

## Means

We will now move over to create Means in R.

We will begin witg a ColMeans in R:

```colMeans(A)
[1] 11 12 13 14 15

```
```colMeans(B)
[1] 27 32 37 42 47

```

Next we will make a RowMeans in R:

```rowMeans(A)
[1]  3  8 13 18 23

```
```rowMeans(B)
[1] 35 36 37 38 39

```

## Sums

Now let us try making Sums in R.

We will begin with ColSums in R:

```colSums(A)
[1] 55 60 65 70 75

```
```colSums(B)
[1] 135 160 185 210 235

```

Now we will try RowSums in R:

```rowSums(A)
[1]  15  40  65  90 115

```
```rowSums(B)

[1] 175 180 185 190 195

```

Lastly we will deal with Advanced matrix operations in R.
First let us try Cholesky Factorization in R, which demands that Matrix should be positive definite:

```S <- matrix(data = c(5, 1, 1, 3), nrow = 2)
R <- chol(S)
t(R) %*% R
[,1] [,2]
[1,]    5    1
[2,]    1    3

```

Now let us deal with Eigenvalues in R:

```y <- eigen(A)
y\$values
[1]  6.864208e+01 -3.642081e+00 -2.749042e-15 -2.127981e-15  8.183097e-16

```
```y\$vectors
[,1]        [,2]        [,3]       [,4]        [,5]
[1,] -0.1079750  0.67495283 -0.03262008  0.1093572 -0.02645084
[2,] -0.2527750  0.36038970 -0.32374283  0.2312142 -0.27585356
[3,] -0.3975750  0.04582657  0.18215364 -0.1969318  0.75630393
[4,] -0.5423751 -0.26873656  0.73740154 -0.7372079 -0.57924382
[5,] -0.6871751 -0.58329969 -0.56319227  0.5935683  0.12524429

```

Let us move over to Moore-Pernrose Generalized Inverse in R:

```library(MASS)
ginv(A)
[,1]   [,2]          [,3]   [,4]   [,5]
[1,] -0.152 -0.096 -4.000000e-02  0.016  0.072
[2,] -0.080 -0.050 -2.000000e-02  0.010  0.040
[3,] -0.008 -0.004 -2.081668e-17  0.004  0.008
[4,]  0.064  0.042  2.000000e-02 -0.002 -0.024
[5,]  0.136  0.088  4.000000e-02 -0.008 -0.056

```
```A %*% ginv(A) %*% A
[,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15
[4,]   16   17   18   19   20
[5,]   21   22   23   24   25

```

Now we will make an QR Decomposition in R:

```y <- qr(A)
y\$qr
[,1]         [,2]          [,3]          [,4]          [,5]
[1,] -29.2403830 -31.12134335 -3.300230e+01 -3.488326e+01 -3.676422e+01
[2,]   0.2051957  -1.20912708 -2.418254e+00 -3.627381e+00 -4.836508e+00
[3,]   0.3761921  -0.03966092 -1.077118e-14 -2.223765e-14 -3.241514e-14
[4,]   0.5471885  -0.43361721  7.421294e-01  2.224166e-15  1.443920e-15
[5,]   0.7181848  -0.82757351  2.061471e-01 -4.638321e-01  1.741866e-15

```
```y\$rank

[1] 2

```
```y\$qraux
[1] 1.034199e+00 1.354295e+00 1.637767e+00 1.885923e+00 1.741866e-15

```
```y\$pivot
[1] 1 2 3 4 5

```

Now let us try making an Inverse in R:

```solve(S)
[,1]        [,2]
[1,]  0.21428571 -0.07142857
[2,] -0.07142857  0.35714286

```
```S %*% solve(S)
[,1]         [,2]
[1,]  1.000000e+00 5.551115e-17
[2,] -2.775558e-17 1.000000e+00

```

Let us try to Solve Linear Equation in R:

```y = Mx
y <- c(1, 2, 3)
M <- matrix(data = c(1, 3, 2, 5, 4, 5, 2, 1, 4), nrow = 3)
x <- solve(M, y)
x
[1]  0.72 -0.20  0.64

```

Now we will make Singular Value Decomposition in R:

```y <- svd(A)
y\$d
[1] 7.425405e+01 3.366820e+00 2.982665e-15 1.539164e-15 2.081370e-16

```
```y\$u
[,1]        [,2]       [,3]        [,4]        [,5]
[1,] -0.09359323 -0.76892152  0.6113052  0.01963160 -0.16099879
[2,] -0.24363203 -0.49055421 -0.5142169 -0.45682855  0.47632830
[3,] -0.39367083 -0.21218690 -0.3594458  0.81870773 -0.02272392
[4,] -0.54370962  0.06618041 -0.1836783 -0.34545623 -0.73954194
[5,] -0.69374842  0.34454772  0.4460359 -0.03605456  0.44693633

```
```y\$v
[,1]       [,2]       [,3]       [,4]        [,5]
[1,] -0.3926228  0.6677180  0.5359956 -0.3354166  0.01429645
[2,] -0.4191312  0.3526033 -0.5060346  0.3600072 -0.56064583
[3,] -0.4456396  0.0374885 -0.5234992 -0.2107002  0.69394096
[4,] -0.4721479 -0.2776263  0.4211198  0.6830453  0.23686977
[5,] -0.4986563 -0.5927410  0.0724184 -0.4969357 -0.38446135

```

Lastly we will make a first order derivative and second derivative hessian matrix:

```library(numDeriv)
sc2.f <- function(x){
n <- length(x)
sum((1:n) * (exp(x) - x)) / n
}
sc2.g <- function(x){
n <- length(x)
(1:n) * (exp(x) - 1) / n
}

x0 <- rnorm(5)
hess <- hessian(func=sc2.f, x=x0)
hessc <- hessian(func=sc2.f, x=x0, "complex")
# Hessian = Jacobian of the gradient
jac <- jacobian(func=sc2.g, x=x0)
jacc <- jacobian(func=sc2.g, x=x0, "complex")
```

Let us display the first order derivative matrix with complex method:

```jacc
[,1]     [,2]      [,3]       [,4]     [,5]
[1,] 0.06255453 0.000000 0.0000000 0.00000000 0.000000
[2,] 0.00000000 2.140093 0.0000000 0.00000000 0.000000
[3,] 0.00000000 0.000000 0.3107356 0.00000000 0.000000
[4,] 0.00000000 0.000000 0.0000000 0.07183351 0.000000
[5,] 0.00000000 0.000000 0.0000000 0.00000000 8.423611

```

Let us also display the second order derivative hessian matrix with Richardson method:

```hess
[,1]         [,2]          [,3]          [,4]          [,5]
[1,]  6.255453e-02 4.218706e-12 -1.377500e-14  5.867048e-12  6.631665e-12
[2,]  4.218706e-12 2.140093e+00  2.577218e-16  2.031748e-12  2.295093e-12
[3,] -1.377500e-14 2.577218e-16  3.107356e-01 -4.314831e-13 -5.864369e-12
[4,]  5.867048e-12 2.031748e-12 -4.314831e-13  7.183351e-02  1.467822e-12
[5,]  6.631665e-12 2.295093e-12 -5.864369e-12  1.467822e-12  8.423611e+00

```

Related Post

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...