Linear Algebra in R

1 Creating Matrix

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

2 Element-wise operations

2.1 Addition

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

2.2 Subtraction

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

2.3 Multiplication

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

2.4 Division

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

2.5 Exponentiation

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

3 Basic Matrix Operations

3.1 Transpose

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

3.2 Multiplication

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

3.3 Outer Product

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
## 
## , , 3, 1
## 
##       [,1]
##  [1,]   23
##  [2,]   46
##  [3,]   69
##  [4,]   92
##  [5,]  115
##  [6,]  138
##  [7,]  161
##  [8,]  184
##  [9,]  207
## [10,]  230
## 
## , , 4, 1
## 
##       [,1]
##  [1,]   24
##  [2,]   48
##  [3,]   72
##  [4,]   96
##  [5,]  120
##  [6,]  144
##  [7,]  168
##  [8,]  192
##  [9,]  216
## [10,]  240
## 
## , , 5, 1
## 
##       [,1]
##  [1,]   25
##  [2,]   50
##  [3,]   75
##  [4,]  100
##  [5,]  125
##  [6,]  150
##  [7,]  175
##  [8,]  200
##  [9,]  225
## [10,]  250
## 
## , , 6, 1
## 
##       [,1]
##  [1,]   26
##  [2,]   52
##  [3,]   78
##  [4,]  104
##  [5,]  130
##  [6,]  156
##  [7,]  182
##  [8,]  208
##  [9,]  234
## [10,]  260
## 
## , , 7, 1
## 
##       [,1]
##  [1,]   27
##  [2,]   54
##  [3,]   81
##  [4,]  108
##  [5,]  135
##  [6,]  162
##  [7,]  189
##  [8,]  216
##  [9,]  243
## [10,]  270
## 
## , , 8, 1
## 
##       [,1]
##  [1,]   28
##  [2,]   56
##  [3,]   84
##  [4,]  112
##  [5,]  140
##  [6,]  168
##  [7,]  196
##  [8,]  224
##  [9,]  252
## [10,]  280
## 
## , , 9, 1
## 
##       [,1]
##  [1,]   29
##  [2,]   58
##  [3,]   87
##  [4,]  116
##  [5,]  145
##  [6,]  174
##  [7,]  203
##  [8,]  232
##  [9,]  261
## [10,]  290
## 
## , , 10, 1
## 
##       [,1]
##  [1,]   30
##  [2,]   60
##  [3,]   90
##  [4,]  120
##  [5,]  150
##  [6,]  180
##  [7,]  210
##  [8,]  240
##  [9,]  270
## [10,]  300

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

3.4 Cross Products

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

3.5 Diagonal

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

4 Combine Matrices

4.1 Horizontally

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

4.2 Vertically

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

5 Means

5.1 ColMeans

colMeans(A)

## [1] 11 12 13 14 15

colMeans(B)

## [1] 27 32 37 42 47

5.2 RowMeans

rowMeans(A)

## [1]  3  8 13 18 23

rowMeans(B)

## [1] 35 36 37 38 39

6 Sums

6.1 ColSums

colSums(A)

## [1] 55 60 65 70 75

colSums(B)

## [1] 135 160 185 210 235

6.2 RowSums

rowSums(A)

## [1]  15  40  65  90 115

rowSums(B)

## [1] 175 180 185 190 195

7 Advanced Matrix Operations

7.1 Cholesky Factorization

Matrix should be positive definite.

https://en.wikipedia.org/wiki/Positive-definite_matrix

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

##      [,1] [,2]
## [1,]    5    1
## [2,]    1    3

7.2 Eigen

https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors

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

7.3 Moore-Pernrose Generalized Inverse

https://en.wikipedia.org/wiki/Moore%E2%80%93Penrose_pseudoinverse

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

7.4 QR Decomposition

https://en.wikipedia.org/wiki/QR_decomposition

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

7.5 Inverse

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

7.6 Solve Linear Equation

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

7.7 Singular Value Decomposition

https://en.wikipedia.org/wiki/Singular_value_decomposition

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

8 Useful Linear Algebra packages

• Matrix package
• matrixStats package
• matlab package