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

Category

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

Advanced Matrix Operations

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

To leave a comment for the author, please follow the link and comment on their blog: R Programming – DataScience+.

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



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers

Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)