-
Notifications
You must be signed in to change notification settings - Fork 15
/
07-matrices-arrays-lists.qmd
929 lines (632 loc) · 29.5 KB
/
07-matrices-arrays-lists.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
---
title: "Matrices, Arrays, and Lists"
author: "Joshua P. French"
date: today
date-format: long
# format: html
format: ipynb
jupyter: ir
execute:
output: false
self-contained: true
title-block-banner: true
bibliography:
- packages_matrices_arrays_lists.bib
---
To open this information in an interactive Colab notebook, click the Open in Colab graphic below.
<a href="https://colab.research.google.com/github/jfrench/DataWrangleViz/blob/master/07-matrices-arrays-lists.ipynb">
<img src="https://colab.research.google.com/assets/colab-badge.svg">
</a>
---
In this module we will cover some additional data structures available in R that were not covered in the Crash Course in R Module. Specifically, we will discuss matrices, arrays, and lists. These data structures are fundamental to performing data analysis in R, particularly if we need to program our own methods.
```{r}
#| include: false
# automatically create a bib database for R packages
knitr::write_bib(c(
.packages(),
'abind',
'microbenchmark'
), 'packages_deeper_dive.bib')
```
::: {.content-hidden when-format="html"}
We start by making sure the necessary R packages are installed.
```{r}
if(!require(microbenchmark, quietly = TRUE)) {
install.packages("microbenchmark",
repos = "https://cran.rstudio.com/")
}
```
```{r}
if(!require(abind, quietly = TRUE)) {
install.packages("abind",
repos = "https://cran.rstudio.com/")
}
```
:::
# Matrices and matrix algebra
## Why are matrices and matrix algebra important knowledge for data scientists? ([YouTube](https://youtu.be/40THrClg9Ho), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=f36ecb67-06e8-45a1-8fb8-b099015ade57))
Understanding how to use matrices is foundational knowledge for any data scientist.
- Matrices are a simple, efficient way of storing many kinds of data.
- Matrix algebra is the most basic way of manipulating data represented in matrix format.
Numeric matrices are often needed to implement complex algorithms used for data science methods related to regression, optimization, and classification.
## What is a matrix? ([YouTube](https://youtu.be/BrNr3cpLxQY), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=35569ba5-8909-4a0b-bb3a-b099015b52e4))
A matrix is a two-dimensional object whose values have the same data type.
- We can have numeric, character, or logical matrices, but only one data type can exist at a time in a matrix.
A matrix $\mathbf{A}$ with $m$ rows and $n$ columns (an $m\times n$ matrix) will be denoted as
$$
\mathbf{A} = \begin{bmatrix}
\mathbf{A}_{1,1} & \mathbf{A}_{1,2} & \cdots & \mathbf{A}_{1,n} \\
\mathbf{A}_{2,1} & \mathbf{A}_{2,2} & \cdots & \mathbf{A}_{2,n} \\
\vdots & \vdots & \ddots & \vdots \\
\mathbf{A}_{m,1} & \mathbf{A}_{m,2} & \cdots & \mathbf{A}_{m,n} \\
\end{bmatrix},
$$
where $\mathbf{A}_{i,j}$ denotes the element in row $i$ and column $j$ of matrix $\mathbf{A}$.
A **column vector** is a matrix with a single column. A **row vector** is a matrix with a single row.
- Vectors are commonly denoted with bold lowercase letters such as $\mathbf{a}$ or $\mathbf{b}$, but this may be simplified to lowercase letters such as $a$ or $b$.
A $p\times 1$ column vector $\mathbf{a}$ may constructed as
$$
\mathbf{a} = [a_1, a_2, \ldots, a_p] =
\begin{bmatrix}
a_1 \\ a_2 \\ \vdots \\ a_p
\end{bmatrix}.
$$
A vector is assumed to be a column vector unless otherwise indicated.
## Creating a matrix ([YouTube](https://youtu.be/OsO39R0S6nY), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=bf3047ad-d8ea-4ba0-99c0-b099015b6066))
A matrix can be created in R by passing a vector to the `matrix` function and specifying the number of rows and columns the matrix will have.
- The values of the vector will be placed into the matrix one column at a time.
- The `nrow` argument specifies the number of rows.
- The `ncol` argument specifies the number of columns.
Let's create a $5\times 2$ matrix with the elements $1, 2, \ldots, 5$ in the first column and $6, 7, \ldots, 10$ in the second column.
```{r}
A <- matrix(1:10, nrow = 5, ncol = 2)
A
```
If we change the `byrow` argument of the `matrix` function to `TRUE`, then the values will be placed row by row, as in the example below.
```{r}
matrix(1:6, nrow = 2, ncol = 3, byrow = TRUE)
```
The `cbind` and `rbind` functions can be used to create a `matrix` object by combining vectors by columns or rows, respectively.
In the example below, we combine the vectors [1, 2, 3, 4] and [5, 6, 7, 8] columnwise to create a matrix.
```{r}
cbind(1:4, 5:8)
```
Alternatively, we can use the `rbind` function to create a matrix by combining the vectors row by row.
```{r}
rbind(1:4, 5:8)
```
The `dim` function can be used to determine/confirm the number of rows or columns in a matrix.
```{r}
dim(A)
```
The `nrow` and `ncol` functions can be used to determine the number of rows and columns, respectively, in a `matrix` object.
```{r}
nrow(A)
```
```{r}
ncol(A)
```
## Subsetting a matrix ([YouTube](https://youtu.be/PEt7HZt7noA), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=4a20bf6d-8b6d-44d9-b304-b09e00c4018a))
We can access a subset of a `matrix` object using the `[` operator and providing vectors with the desired rows and columns separated by a comma.
If `A` is a matrix, then we can select the subset of `A` using the syntax `A[r, c]`, where `r` and `c` are `numeric` vector indicating the desired rows and columns, respectively, that we want to select.
- If we select a single row or column, then the subsetted object will simplify to `vector` object instead of a `matrix`.
- If we don't supply a vector of rows or columns, then all rows or columns will be returned.
- We can use the `-` notation to select all but the specified rows or columns. This is similar to the syntax used for data frames.
This is easiest to understand by walking through a set of examples.
Let's define a $4\times 3$ matrix, `B`, below.
```{r}
B <- matrix(1:12, ncol = 3)
B
```
Let's select columns 2 and 3 of `B`.
```{r}
B[, 2:3]
```
Let's select rows 1 and 4 of `B`.
```{r}
B[c(1, 4), ]
```
Let's select columns 2 and 3 of rows 1 and 4 of `B`.
```{r}
B[c(1, 4), 2:3]
```
Let's select all columns of `B` except the third.
```{r}
B[,-3]
```
Our matrix simplifies to a vector if we only select a single row or column.
We select the second row of `B`.
```{r}
B[2, ]
```
We select the second column of `B`.
```{r}
B[,2]
```
If we want to retain the `matrix` structure of `B` when we select a single row or column, then we set the optional `drop` argument to `FALSE` after the column argument.
```{r}
B[2, , drop = FALSE]
```
```{r}
B[, 2, drop = FALSE]
```
If we want to extract specific elements of a `matrix` object and not a rectangular subset of the object, then we can use a `matrix` indicating the elements we want to extract inside `[`.
Let's say we want to extract the elements in positions $(1, 2)$, $(3, 1)$, and $(4, 3)$. We create a matrix describing those positions. Each row of the matrix describes an element we want to extract.
```{r}
pos <- cbind(c(1, 3, 4), c(2, 1, 3))
pos
```
We use `pos` to extract those elements from `B`.
```{r}
B[pos]
```
## Matrix operations
### Addition and subtraction ([YouTube](https://youtu.be/ChiorfijWqY), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=f0afa263-b765-4054-b808-b09e0159ccb8))
We can add or subtract the elements of two matrices as long as they have the same size.
We create two $2\times 3$ matrices below.
```{r}
A <- matrix(1:6, nrow = 2, ncol = 3, byrow = TRUE)
A
```
```{r}
B <- matrix(c(2, 1, 9, 3, 1, 1), nrow = 2, ncol = 3)
B
```
We can add the two matrices below using the `+` operator. The result is simply the sum of the values in the same position of the matrices.
```{r}
A + B
```
Subtraction works in the same way, as shown below.
```{r}
A - B
```
### Scalar matrix multiplication ([YouTube](https://youtu.be/gCB4EWMuF0Y), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=ec1b70ff-5a35-4c99-81b6-b09e015a8b6b))
A matrix multiplied by a scalar value $c\in\mathbb{R}$ is the matrix obtained by multiplying each element of the matrix by $c$.
To perform scalar multiplication in R, we simply use the `*` operator to multiply the scalar value by the matrix.
In the code below, we double the values in `A`.
```{r}
2 * A
```
### Matrix multiplication ([YouTube](https://youtu.be/QRTMoCoqai4), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=377475fa-ea2c-4b8c-8daf-b09e015b5d8e))
Consider two matrices $\mathbf{A}$ and $\mathbf{B}$. The matrix product $\mathbf{AB}$ is only defined if the number of columns in $\mathbf{A}$ matches the number of rows in $\mathbf{B}$.
Assume $\mathbf{A}$ is an $m\times n$ matrix and $\mathbf{B}$ is an $n\times p$ matrix. $\mathbf{AB}$ will be an $m\times p$ matrix and
$$
(\mathbf{AB})_{i,j} = \sum_{k=1}^{n} \mathbf{A}_{i,k}\mathbf{B}_{k,j}.
$$
We consider multiplying the two matrices below.
$$
\begin{bmatrix}
1 & 2 & 3 \\
4 & 5 & 6
\end{bmatrix}
\begin{bmatrix}
1 & 4\\
2 & 5\\
3 & 6
\end{bmatrix}=
\begin{bmatrix}
1\cdot 1 + 2 \cdot 2 + 3 \cdot 3 & 1 \cdot 4 + 2 \cdot 5 + 3 \cdot 6\\
4\cdot 1 + 5 \cdot 2 + 6 \cdot 3 & 4 \cdot 4 + 5 \cdot 5 + 6 \cdot 6\\
\end{bmatrix}=
\begin{bmatrix}
14 & 32\\
32 & 77\\
\end{bmatrix}.
$$
We multiply two matrices of the correct size using the `%*%` operator. Let's multiply these same matrices using R
```{r}
A
```
```{r}
B <- matrix(1:6, nrow = 3, ncol = 2)
B
```
```{r}
A %*% B
```
### Elementwise matrix multiplication ([YouTube](https://youtu.be/DSqjWcorneI), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=25239575-20f2-4e12-a35a-b09e015c3588))
Consider two matrices $\mathbf{A}$ and $\mathbf{B}$ that have the same dimensions. The elementwise product, denoted $\mathbf{A} \circ \mathbf{B}$ is the product of the elements in the same position of the matrix.
Assume $\mathbf{A}$ and $\mathbf{B}$ are $n\times m$ matrices. $\mathbf{A} \circ \mathbf{B}$ will be an $n\times m$ matrix and
$$
(\mathbf{A} \circ \mathbf{B})_{i,j} = \mathbf{A}_{i,j}\mathbf{B}_{i,j}.
$$
We perform elementwise multiplication of the two matrices below.
$$
\begin{bmatrix}
1 & 4\\
2 & 5\\
3 & 6
\end{bmatrix} \circ
\begin{bmatrix}
7 & 10\\
8 & 11\\
9 & 12
\end{bmatrix}=
\begin{bmatrix}
1 \cdot 7 & 4 \cdot 10\\
2 \cdot 8 & 5 \cdot 11\\
3 \cdot 9 & 6 \cdot 12
\end{bmatrix}=
\begin{bmatrix}
7 & 40\\
16 & 55\\
27 & 72
\end{bmatrix}.
$$
We perform elementwise multiplication of two matrices of the correct size using the `*` operator. Let's multiply these same matrices using R
```{r}
A <- matrix(1:6, ncol = 2)
A
```
```{r}
B <- matrix(7:12, ncol = 2)
```
```{r}
B
```
```{r}
A * B
```
### Transpose ([YouTube](https://youtu.be/R_v94nP4O70), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=0fb7b469-0f72-407d-bf10-b09e015d2526))
The **transpose** of a matrix $\mathbf{A}$, denoted $\mathbf{A}^T$, exchanges the rows and columns of the matrix. More formally, the $i,j$ element of $\mathbf{A}^T$ is the $j,i$ element of $\mathbf{A}$, i.e., $(\mathbf{A}^T)_{i,j} = \mathbf{A}_{j,i}$.
We can use the `t` function to obtain the transpose of a matrix.
We print the matrix `A` that we have previously defined.
```{r}
A
```
We compute the transpose of `A` using `t`.
```{r}
t(A)
```
## Special matrices ([YouTube](https://youtu.be/-05Z3BG3ZZA), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=5d546ef4-c64f-471d-b1a4-b09e015e3620))
### Square matrices
A matrix is **square** if the number of rows equals the number of columns.
The **diagonal elements** of an $n\times n$ square matrix $\mathbf{A}$ are the elements $\mathbf{A}_{i,i}$ for $i = 1, 2, \ldots, n$.
Non-diagonal elements of $\mathbf{A}$ are called **off-diagonal** elements.
### Identity matrix
The $n\times n$ **identity** matrix, $\mathbf{I}_{n \times n}$, is a matrix for which the diagonal elements are all 1 and the non-diagonal elements are all 0.
- Context often makes it clear what the dimensions of an identity matrix are, so the notation for $\mathbf{I}_{n\times n}$ is often simplified to $\mathbf{I}$ or $I$.
The syntax `diag(n)` will create an $n\times n$ identity matrix in R.
We create a $3\times 3$ identify matrix below.
```{r}
diag(3)
```
### Inverse matrix
An $n\times n$ matrix $\mathbf{A}$ is invertible if there exists a matrix $\mathbf{B}$ such that $\mathbf{AB}=\mathbf{BA}=\mathbf{I}_{n\times n}$.
The inverse of $\mathbf{A}$ is denoted $\mathbf{A}^{-1}$.
- Inverse matrices only exist for square matrices.
The `solve` function can be used to compute the inverse of matrix, e.g., `solve(A)`.
- This will likely cause numeric issues, so don't do this.
Inverse matrices are nearly always multiplied by a different matrix, e.g., $\mathbf{b}$.
To compute $\mathbf{A}^{-1} \mathbf{b}$, we use the syntax `solve(A, b)`.
In the example below, we use `solve` to compute $\mathbf{A}^{-1} \mathbf{b}$.
```{r}
A <- matrix(c(3, 0.5, 0.2,
0.5, 2, 0.1,
0.2, 0.1, 1),
nrow = 3)
A
```
```{r}
b <- c(1, 0.7, 3)
b
```
```{r}
solve(A, b)
```
Because vectors are assumed to be column vectors, the (rounded) result above is equivalent to
$$
\begin{bmatrix}
0.11\\
0.18\\
2.96
\end{bmatrix}.
$$
## Example: Computing the line of best fit using matrix algebra ([YouTube](https://youtu.be/rSVxdqRCkuE), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=4d6cbd52-1169-49fd-b9a1-b09e015eb369))
A common application of matrix algebra is finding the "line of best fit" for the points observed in a scatter plot.
A line can be written as
$$y = a + b x,$$
where:
- $a$ is the y-intercept (the value of $y$ when $x = 0$).
- $b$ is the slope of the line. $y$ changes by $b$ when $x$ increases by 1.
Suppose that we observe $n$ paired sets of $(x, y)$ values
$$(x_1, y_1), (x_2, y_2), \ldots, (x_n, y_n).$$
We want to find the values of $a$ and $b$ that minimize the residual sum of squares function
$$
RSS(a, b) = \sum_{i=1}^n (y_i - a - b x_i)^2.
$$
We can find the values of $a$ and $b$ using matrix algebra. First, we must define a matrix and a column vector.
Define
$$
\mathbf{X} =
\begin{bmatrix}
1 & x_1 \\
1 & x_2 \\
\vdots & \vdots \\
1 & x_n
\end{bmatrix}.
$$
Define
$$
\mathbf{y} =
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots \\
y_n
\end{bmatrix}.
$$
The values of $a$ and $b$ that minimize $RSS(a, b)$, which we will call $\hat{a}$ and $\hat{b}$, are computed through the equation
$$
\begin{bmatrix}
\hat{a} \\
\hat {b}
\end{bmatrix} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{y}.
$$ {#eq-ols-solution}
Let says we have the following 6 pairs of $(x, y)$ values:
$$(1, 4.9), (1.3, 4.2), (1.7, 4.0), (2.1, 5.7), (3, 4.6), (3.5, 6.8).$$
We turn these into two vectors, `x` and `y`, below.
```{r}
x <- c(1, 1.3, 1.7, 2.1, 3, 3.5)
y <- c(4.9, 4.2, 4.0, 5.7, 4.6, 6.8)
```
We create a scatter plot of the vectors below.
```{r}
plot(x, y)
```
We now create the $\mathbf{X}$ matrix using the `cbind` function by combining a vector of 1s with `x`.
```{r}
X <- cbind(1, x)
X
```
Now, we use @eq-ols-solution to determine the line of best fit.
For clarity, we'll first compute $\mathbf{X}^T \mathbf{X}$, then $\mathbf{X}^T \mathbf{y}$, and then use the `solve` function to obtain our solution.
```{r}
xtx <- t(X) %*% X
xty <- t(X) %*% y
solve(xtx, xty)
```
Thus, our solution for $\hat{a}$ and $\hat{b}$ that minimizes $RSS(a, b)$ is $\hat{a} = 3.58$ and $\hat{b} = 0.69$.
We use the `abline` function to add the line of best fit to our scatter plot.
```{r}
plot(x, y)
abline(a = 3.58, b = 0.69)
```
## Tips for optimizing matrix algebra operations ([YouTube](https://youtu.be/cBovzIQvw78), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=5e2bb274-0855-40ef-a60e-b09e015fec12))
It is straightforward to perform matrix algebra in R. However, it is more difficult to write *efficient* (faster) matrix code. We provide some basic tips below.
### `crossprod` and `tcrossprod`
The `crossprod` and `tcrossprod` functions are faster approaches to performing matrix multiplication when one of the matrices must first be transposed.
If we have matrices `A` and `B` of compatible dimensions, then:
- `t(A) %*% B` is equivalent to `crossprod(A, B)`.
- `A %*% t(B)` is equivalent to `tcrossprod(A, B)`.
`crossprod` and `tcrossprod` are faster than their equivalent approaches because the equivalent approaches compute the transpose of matrix in memory before performing the multiplication. The `crossprod` and `tcrossprod` functions modify their multiplication algorithms so that the transpose matrix doesn't have to be created prior to multiplication.
To illustrate this difference, we create `A`, a $50\times 2$ matrix of values and `B`, a $50 \times 2$ matrix of values by drawing the appropriate number of values from a standard normal distribution. We use the `set.seed` function to make the example reproducible.
```{r}
set.seed(314)
A <- matrix(rnorm(100), ncol = 2)
B <- matrix(rnorm(100), ncol = 2)
```
First, we use the `all.equal` function to confirm that the standard and optimized approaches produce the same results. `all.equal` confirms that the supplied objects are equal to each other.
```{r}
all.equal(t(A) %*% B, crossprod(A, B))
```
```{r}
all.equal(A %*% t(B), tcrossprod(A, B))
```
We will use the `microbenchmark` function from the **microbenchmark** package [@R-microbenchmark] to time the optimized and standard approaches to computing $\mathbf{A}^T \mathbf{B}$ and $\mathbf{A} \mathbf{B}^T$.
We start by loading the **microbenchmark** package.
```{r}
library(microbenchmark)
```
We now compute $\mathbf{A}^T \mathbf{B}$ using the standard approach and using `crossprod` within the `microbenchmark` function. The `microbenchmark` function will perform each computation 100 times by default to determine the typical speed. We then use the `plot` function to plot the results.
```{r}
t1 <- microbenchmark::microbenchmark(t(A) %*% B, crossprod(A, B))
plot(t1)
```
While the results shown in the above plot will depend on the speed of the computer used to perform the calculation, the boxplot shown for `crossprod(A, B)` should be lower then `t(A) %*% B`, which means the `crossprod` approach is faster than the standard approach.
We perform a simialr analysis for computing $\mathbf{A}\mathbf{B}^T$.
```{r}
t2 <- microbenchmark(A %*% t(B), tcrossprod(A, B))
plot(t2)
```
Once again, the plot of our timing results should show that `tcrossprod(A, B)` can be computed faster than `A %*% t(B)`.
### Order of matrix multiplication
Matrix multiplication is associative, meaning that if we have compatible matrices $\mathbf{A}$, $\mathbf{B}$ and $\mathbf{C}$, then $$(\mathbf{A B})\mathbf{C} = \mathbf{A} (\mathbf{B C}).$$
The order of matrix multiplication can greatly impact how long it takes to perform the operations.
**Prioritize performing operations that produce smaller matrices first.** Smaller can mean the total number of rows/columns or long/thin matrices, depending on the context.
We create three matrices below filled with random values: `A` is a matrix of size $500\times 2$, `B` is a matrix of size $2 \times 50$, `C` is a $50\times 1$ column vector.
```{r}
set.seed(1999)
A <- matrix(rnorm(1000), ncol = 2)
B <- matrix(rnorm(100), nrow = 2)
C <- rnorm(50)
```
Computing $\mathbf{A}(\mathbf{BC})$ is much faster than computing $(\mathbf{AB})\mathbf{C}$ because $\mathbf{BC}$ first produces $2\times 1$ matrix while $\mathbf{AB}$ produces a much larger $500 \times 50$ matrix.
We confirm this by timing these operations using the `microbenchmark` function and plotting the results.
```{r}
t3 <- microbenchmark((A %*% B) %*% C,
A %*% (B %*% C))
plot(t3, las = 2, cex.axis = 0.5)
```
# Arrays
## What is an array? ([YouTube](https://youtu.be/Gd5Z5TR9rQI), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=1d53c909-ca7a-4ae6-b63a-b09e0161a2d3))
An array is a multidimensional, homogeneous (one data type) data structure.
- Arrays can be 1, 2, 3, $\ldots$ dimensions.
- Arrays are a generalization of a matrix.
We will discuss 3-dimensional arrays because they are the easiest to visualize.
A 3-dimensional array is like a loaf of bread.
- Each slice of bread is a two-dimensional matrix/array.
- We stack the matrix slices side-by-side to get the third dimension of our array.
::: {.content-hidden when-format="ipynb"}
![[Sliced bread](https://www.pexels.com/photo/wheat-bread-slices-166021/)](./figures/sliced-bread.jpg)
:::
::: {.content-hidden when-format="html"}
![Sliced bread](https://images.pexels.com/photos/166021/pexels-photo-166021.jpeg)
:::
## Creating an array ([YouTube](https://youtu.be/rI-S4DgESLU), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=78e4896e-a7bc-42fb-99bf-b09e0162500c))
We can create an array in R using the `array` function.
The two main arguments of the `array` function are:
- `data`: an atomic vector.
- `dim`: a numeric vector indicating the size of each dimension.
Similar to the `matrix` function, the `array` function fills in the columns of each matrix before repeating the process for each element of the third dimension.
We create a 3-dimensional array with dimension $3\times 2 \times 4$ below.
```{r}
A <- array(1:24, dim = c(3, 2, 4))
print(A)
```
Alternatively, we can bind matrices or arrays together using the `abind` function in the **abind** package [@R-abind].
This is actually a very complex function, but we illustrate its use for creating a 3-dimensional array by binding 2-dimensional matrices.
- Run `?abind::abind` in the Console for more details.
We first load the **abind** package.
```{r}
library(abind)
```
We now bind two matrices along a 3rd-dimension to create a 3-dimensional array.
```{r}
A <- abind(
matrix(1:4, nrow = 2),
matrix(5:8, nrow = 2),
along = 3)
print(A)
```
We create a second 3-dimensional array (of size $2\times 2 \times 2$ named `B` below.
```{r}
B <- array(9:16, dim = c(2, 2, 2))
```
To bind `A` and `B` together, we can use the `abind` function.
```{r}
AB <- abind(A, B)
print(AB)
```
The `dim` function can be used to determine the size of an array.
```{r}
dim(A)
```
## Subsetting an array ([YouTube](https://youtu.be/tfynOZxThHA), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=dc55a667-f106-4254-a1ed-b09e0162f98f))
We can subset an array in much the same way as a matrix. We use the `[` operator and provide vectors indicating the desired positions in each dimension that we want to extract, with each dimension being separated by a comma.
If `A` is an array with 3 dimensions, then we can select a rectangular subset of `A` using the syntax `A[v1, v2, v3]`, where `v1`, `v2`, and `v3` are `numeric` vectors indicating the desired elements in each dimension, respectively, that we want to select.
We provide some examples below.
```{r}
A <- array(1:24, dim = c(3, 2, 4))
print(A)
```
Here, we subset the first element of the first dimension of `A`. This will run across the two elements of the second dimension of `A` and across all 4 elements of the third dimension.
```{r}
print(A[1, , ])
```
We subset the second element of the second dimension of `A`. This will run across all 3 elements of the first dimension of `A` and across all 4 elements of the 3rd dimension.
```{r}
print(A[, 2, ])
```
We subset the second and third elements of the third dimension of `A`, i.e., the second and third slices of the third dimension of `A`. This will run across all 3 elements of the first dimension of `A` and across both elements of the second dimension.
```{r}
print(A[, , 2:3])
```
If we subset an array such that the subset can be thought of as a vector (e.g., a $1\times 1 \times 4$ array is like a vector with length 4), then R will automatically convert the subsetted object to a vector because the `[` operator has a hidden `drop` argument that by default coerces the result to the lowest possible dimension.
```{r}
print(A[1, 1, ])
```
To avoid this behavior, we must set `drop` to `FALSE`.
```{r}
print(A[1, 1, , drop = FALSE])
```
If we want to select specific *elements* of an array (and not a rectangular subset) then we can pass a matrix indicating the elements we want to select.
- The number of columns of the matrix needs to match the number of dimensions of the array we are subsetting.
- Each row of the matrix indicates an element we want to select.
We create an matrix, `v`, to subset the elements in positions $(1, 1, 1)$ and $(3, 2, 4)$ of `A`.
```{r}
v <- cbind(c(1, 3), c(1, 2), c(1, 4))
v
```
We now use `v` to subset these elements of `A`.
```{r}
print(A[v])
```
We can use the `-` syntax we learned with `matrix` objects to subset parts of an array that are not in certain positions, but we don't demonstrate its usage here.
# Lists
## What is a list? ([YouTube](https://youtu.be/VGqZg0Na7kU), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=ce680bde-c6d6-480d-a331-b09e01643cb4))
A list in R is a non-atomic vector.
- It is a vector, meaning it is one-dimensional.
- It is non-atomic, meaning it can hold data beyond the basic types (`character`, `integer`, `double`, `logical`, `complex`, `raw`).
Lists are commonly used in two different ways:
1. A vector with elements of the same non-atomic type, e.g., a vector of `lm` objects from the `lm` function.
2. A vector of elements of many different types. This is often used in objected-oriented programming in R, which we will not discuss.
We will discuss the first type of list in the context of loops and apply functions, where they make more sense.
We will discuss the second type below.
## Creating a list ([YouTube](https://youtu.be/l5-qDCntYZg), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=b83d4606-7365-4576-9e5a-b09e017310d9))
A list is stored as a `list` object.
A `list` object can be created using the `list` function. Simply pass the objects you want to include the `list` as arguments to the `list` function.
We provide an example below. We create a list that includes a vector, a matrix, and a function.
```{r}
list(1:10, matrix(1:10, nrow = 2), mean)
```
It is common to assign the elements of a list names to make them easier to reference. To do that, we simply use the syntax `name = x` in our list, where `name` is the name we want to give the element and `x` is the object we are including in our list.
```{r}
mylist <- list(a = 1:10, b = matrix(1:10, nrow = 2), f = mean)
mylist
```
We can create a new list by combining two or more other lists using the `c` function.
We combine 3 simple lists below using the `c` function.
```{r}
c(list(1, 1:2), list("a", c("a", "b")), list(TRUE, c(TRUE, FALSE)))
```
## Subsetting a list ([YouTube](https://youtu.be/mc5oUxvSVZM), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=fdd25b1c-f02e-4ba2-8231-b09e0173c626))
Subsetting a list is more similar to subsetting an atomic vector than it is to subsetting a matrix or array.
We create a simple 4-element list with named` components below.
```{r}
a <- list(a = 1, b = 1:2, c = 1:3, d = 1:4)
```
To subset elements of the list, we append `[idx]` to the end of the list's name, where `idx` is a vector with the indices of the elements we want to extract or a logical vector indicating the elements we want to extract.
In the example below, we extract the 1st and 3rd elements of `a`.
```{r}
a[c(1, 3)]
```
In the next example, we use a logical vector to subset the 2nd and 4th elements of `a`.
```{r}
a[c(FALSE, TRUE, FALSE, TRUE)]
```
Since the elements of the list are named, we can also use a vector of character strings to indicate the elements we want to subset.
We subset elements `c` and `d` below.
```{r}
a[c("c", "d")]
```
If we want to extract a specific element out of our list, then we append `[[i]]` to the list's name, with `i` indicating the index of the element we want to extract.
In the example below, we extract the second element of `a`. Notice that the result is the vector `1:2`, which is what is contained in the 2nd element of `a`.
```{r}
a[[2]]
```
What is the difference between the previous result and what we get from the following command?
```{r}
a[2]
```
In the former example, we extracted the object that was in the 2nd element of `a`. In the latter example, we get a list with only one element. The only element of this list is what was previously the 2nd element of `a`.
For a list with named components, we can use the `$` operator to extract a specific element.
In the code below, we extract the object in element `"c"` of `a`.
```{r}
a$c
```
## Nested lists ([YouTube](https://youtu.be/wkgPuUgtPR8), [Panopto](https://ucdenver.hosted.panopto.com/Panopto/Pages/Viewer.aspx?id=b74d8590-6e6d-4ef7-a541-b09e017450b9))
We can create lists of lists of lists, which created "nested" lists because the lists are nested in each other.
In that case, to subset elements of the lists, we might need to stack usage of `[]` or `[[]]` side by side to subset or extract the desired elements.
We create a nested list below.
```{r}
b <- list(list(1, 1:2, list(3, 3:4)), 5:7)
b
```
The first element of the top-level of the list is another list with the second element is the vector `5:7`.
```{r}
b[[1]]
```
```{r}
b[[2]]
```
`b[[1]]` is a list with 3 elements: the first element is the vector `1`, the second element is the vector `[3]`, while the third element is another list with two elements.
We can subset the first two elements of `b[[1]]` by appending `[1:2]` to the previous syntax.
```{r}
b[[1]][1:2]
```
We can extract the third element of `b[[1]]` by appending `[[3]]` to the previous syntax.
```{r}
b[[1]][[3]]
```
We can access information contained in the 3rd element of the list contained in element of `b` in a similar way.